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Abstract. We study the general problem of minimizing a convex function over a compact 
O convex domain. We will investigate a simple iterative approximation algorithm based on the 
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function over a class of diagonally dominant symmetric matrices. 

We show that our proposed first-order method also applies to nuclear norm and max-norm 
I matrix optimization problems. For nuclear norm regularized optimization, such as matrix com- 

QQ pletion and low-rank recovery, we demonstrate the practical efficiency and scalability of our 

algorithm for large matrix problems, as e.g. the Netflix dataset. For general convex optimiza- 
tion over bounded matrix max-norm, our algorithm is the first with a convergence guarantee, 
Kjjl to the best of our knowledge. 

• ^ 

(This article consists of the first two chapters of the author's PhD thesis [.Jagll].) 



method by Frank & Wolfe [FW56], that does not need projection steps in order to stay inside 
the optimization domain. Instead of a projection step, the linearized problem defined by a 
current subgradient is solved, which gives a step direction that will naturally stay in the do- 
main. Our framework generalizes the sparse greedy algorithm of [FW56] and its primal-dual 
analysis by [ClalO] (and the low-rank SDP approach by [Haz08]) to arbitrary convex domains. 
Analogously, we give a convergence proof guaranteeing e-small duality gap after O(^) iterations. 

The method allows us to understand the sparsity of approximate solutions for any ^i-regularized 
convex optimization problem (and for optimization over the simplex), expressed as a function of 
the approximation quality. We obtain matching upper and lower bounds of 0(^) for the sparsity 
for ^i-problems. The same bounds apply to low-rank semidefinite optimization with bounded 
trace, showing that rank O(^) is best possible here as well. As another application, we obtain 
sparse matrices of O(-) non-zero entries as e-approximate solutions when optimizing any convex 
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1 Introduction 



Motivation. For the performance of large scale approximation algorithms for convex optimiza- 
tion, the trade-off between the number of iterations on one hand, and the computational cost 
per iteration on the other hand, is of crucial importance. The lower complexity per iteration 
is among the main reasons why first-order methods (i.e., methods using only information from 
the first derivative of the objective function), as for example stochastic gradient descent, are 
currently used much more widely and successfully in many machine learning applications — de- 
spite the fact that they often need a larger number of iterations than for example second-order 
methods. 

Classical gradient descent optimization techniques usually require a projection step in each 
iteration, in order to get back to the feasible region. For a variety of applications, this is a 
non-trivial and costly step. One prominent example is semidefinite optimization, where the 
projection of an arbitrary symmetric matrix back to the PSD matrices requires the computation 
of a complete eigenvalue-decomposition. 

Here we study a simple first-order approach that does not need any projection steps, and is 
applicable to any convex optimization problem over a compact convex domain. The algorithm 
is a generalization of an existing method originally proposed by Frank & Wolfe [FW56], which 
was recently extended and analyzed in the seminal paper of Clarkson [ClalO] for optimization 
over the unit simplex. 

Instead of a projection, the primitive operation of the optimizer here is to minimize a linear ap- 
proximation to the function over the same (compact) optimization domain. Any (approximate) 
minimizer of this simpler linearized problem is then chosen as the next step-direction. Because 
all such candidates are always feasible for the original problem, the algorithm will automatically 
stay in our convex feasible region. The analysis will show that the number of steps needed is 
roughly identical to classical gradient descent schemes, meaning that O (^) steps suffice in order 
to obtain an approximation quality of e > 0. 

The main question about the efficiency per iteration of our algorithm, compared to a classical 
gradient descent step, can not be answered generally in favor of one or the other. Whether a 
projection or a linearized problem is computationally cheaper will crucially depend on the shape 
and the representation of the feasible region. Interestingly, if we consider convex optimization 
over the Euclidean ||.||2-ball, the two approaches fully coincide, i.e., we exactly recover classical 
gradient descent. However there are several classes of optimization problems where the lin- 
earization approach we present here is definitely very attractive, and leads to faster and simpler 
algorithms. This includes for example ^i-regularized problems, which we discuss in Sections 4 
and 5, as well as semidefinite optimization under bounded trace, as studied by [Haz08], see 
Section 7. 

Sparsity and Low-Rank. For these mentioned specific classes of convex optimization problems, 
we will additionally demonstrate that our studied algorithm leads to (optimally) sparse or low- 
rank solutions. This property is a crucial side-effect that can usually not be achieved by classical 
optimization techniques, and corresponds to the coreset concept known from computational 
geometry, see also [GJ09]. More precisely, we show matching upper and lower bounds of (^) 
for the sparsity of solutions to general £i-regularized problems, and also for optimizing over 
the simplex, if the required approximation quality is e. For matrix optimization, an analogous 
statement will hold for the rank in case of nuclear norm regularized problems. 

Applications. Applications of the first mentioned class of ^i-regularized problems do include 
many machine learning algorithms ranging from support vector machines (SVMs) to boosting 
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and multiple kernel learning, as well as £2-support vector regression (SVR), mean- variance anal- 
ysis in portfolio selection [Mar52], the smallest enclosing ball problem [BC07], £i -regularized 
least squares (also known as basis pursuit de- noising in compressed sensing), the Lasso [Tib96], 
and £i-regularized logistic regression [KKB07] as well as walking of artificial dogs over rough 
terrain [KBP+10]. 

The second mentioned class of matrix problems, that is, optimizing over semidefinite matrices 
with bounded trace, has applications in low-rank recovery [FHBOl, CR09, CTIO], dimension- 
ality reduction, matrix factorization and completion problems, as well as general semidefinite 
programs (SDPs). 

Further applications to nuclear norm and max-norm optimization, such as sparse/robust PCA 
will be discussed in Section 11. 

History and Related Work. The class of first-order optimization methods in the spirit of Frank 
and Wolfe [FW56] has a rich history in the literature. Although the focus of the original paper 
was on quadratic programming, its last section [FW56, Section 6] already introduces the general 
algorithm for minimizing convex functions using the above mentioned linearization idea, when 
the optimization domain is given by linear inequality constraints. In this case, each intermediate 
step consists of solving a linear program. The given convergence guarantee bounds the primal 
error, and assumes that all internal problems are solved exactly. 

Later [DH78] has generalized the same method to arbitrary convex domains, and improved 
the analysis to also work when the internal subproblems are only solved approximately, see 
also [DunSO]. Patriksson in [Pat93, Pat98] then revisited the general optimization paradigm, 
investigated several interesting classes of convex domains, and coined the term "cost approxi- 
mation" for this type of algorithms. More recently, [Zha03] considered optimization over convex 
hulls, and studies the crucial concept of sparsity of the resulting approximate solutions. However, 
this proposed algorithm does not use linear subproblems. 

The most recent work of Clarkson [ClalO] provides a good overview of the existing lines of re- 
search, and investigates the sparsity solutions when the optimization domain is the unit simplex, 
and establishing the connection to coreset methods from computational geometry. Furthermore, 
[ClalO] was the first to introduce the stronger notion of convergence in primal-dual error for this 
class of problems, and relating this notion of duality gap to Wolfe duality. 

Our Contributions. The character of this article mostly lies in reviewing, re-interpreting and 
generalizing the existing approach given by [ClalO], [Haz08] and the earlier papers by [Zha03, 
DH78, FW56], who do deserve credit for the analysis techniques. Our contribution here is to 
transfer these methods to the more general case of convex optimization over arbitrary bounded 
convex subsets of a vector space, while providing stronger primal-dual convergence guarantees. 
To do so, we propose a very simple alternative concept of optimization duality, which will allow 
us to generalize the stronger primal-dual convergence analysis which [ClalO] has provided for 
the the simplex case, to optimization over arbitrary convex domains. So far, no such guarantees 
on the duality gap were known in the literature for the Frank- Wolfe- type algorithms [FW56], 
except when optimizing over the simplex. Furthermore, we generalize Clarkson's analysis [ClalO] 
to work when only approximate linear internal optimizers are used, and to arbitrary starting 
points. Also, we study the sparsity of solutions in more detail, obtaining upper and lower bounds 
for the sparsity of approximate solutions for a wider class of domains. 

Our proposed notion of duality gives simple certificates for the current approximation quality, 
which can be used for any optimization algorithm for convex optimization problems over bounded 
domain, even in the case of non-differentiable objective functions. 
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We demonstrate the broad applicability of our general technique to several important classes 
of optimization problems, such as ii- and ^oo-regularized problems, as well as semidefinite opti- 
mization with uniformly bounded diagonal, and sparse semidefinite optimization. 

Later in Section 11 we will give a simple transformation in order to apply the first-order 
optimization techniques we review here also to nuclear norm and max-norm matrix optimization 
problems. 

Acknowledgments. Credit for the important geometric interpretation of the duality gap over 
the spectahedron as the distance to the linearization goes to Soeren Laue. Furthermore, the 
author would like to thank Marek Sulovsky, Bernd Gartner, Arkadi Nemirovski, Elad Hazan, 
Joachim Giesen, Sebastian Stich, Michel Baes, Michael Biirgisser and Christian Lorenz Miiller 
for helpful discussions and comments. 

2 The Poor Man's Approach to Convex Optimization and Duality 

The Idea of a Duality given by Supporting Hyperplanes. Suppose we are given the task of 
minimizing a convex function / over a bounded convex set D C M", and let us assume for the 
moment that / is continuously differentiable. 

Then for any point x £ D, it seems natural to consider the tangential "supporting" hyperplane 
to the graph of the function / at the point {x, f{x)). Since the function / is convex, any such 
linear approximation must lie below the graph of the function. 

Using this linear approximation for each point x £ D, we define a dual function value uj{x) 
as the minimum of the linear approximation to / at point x, where the minimum is taken over 
the domain D. We note that the point attaining this linear minimum also seems to be good 
direction of improvement for our original minimization problem given by /, as seen from the 
current point x. This idea will lead to the optimization algorithm that we will discuss below. 

As the entire graph of / lies above any such linear approximation, it is easy to see that 
oj{x) < f{y) holds for each pair x,y £ D. This fact is called weak duality in the optimization 
literature. 

This rather simple definition already completes the duality concept that we will need in this 
paper. We will provide a slightly more formal and concise definition in the next subsection, 
which is useful also for the case of non-differentiable convex functions. The reason we call this 
concept a poor man's duality is that we think it is considerably more direct and intuitive for the 
setting here, when compared to classical Lagrange duality or Wolfe duality, see e.g. [BV04]. 

2.1 Subgradients of a Convex Function 

In the following, we will work over a general vector space X equipped with an inner product 
(.,.). As the most prominent example in our investigations, the reader might always think of 
the case = M" with {x,y) = x'^y being the standard Euclidean scalar product. 

We consider general convex optimization problems given by a convex function / : — >• M over 
a compact^ convex domain D C X, or formally 



In order to develop both our algorithm and the notion of duality for such convex optimization 
problems in the following, we need to formally define the supporting hyperplanes at a given 

^Here we call a set D Q X compact if it is closed and bounded. See [KZ05] for more details. 
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point X £ D. These planes coincide exactly with the well-studied concept of subgradients of a 
convex function. 

For each point x £ D, the subdifferential at x is defined as the set of normal vectors of the 
affine linear functions through {x,f{x)) that lie below the function /. Formally 

df{x) := {d^eX \ f{y) > fix) + {y-x, 4) Vy G 7^} . (2) 

Any element dx G df{x) is called a subgradient to / at x. Note that for each x, df{x) is a 
closed convex set. Furthermore, if / is differentiable, then the subdifferential consists of exactly 
one element for each x G D, namely df{x) = {V/(x)}, as explained e.g. in [Ncm05, KZ05]. 

If we assume that / is convex and lower semicontinuous^ on D, then it is known that df{x) 
is non-empty, meaning that there exists at least one subgradient dx for every point x £ D. 
For a more detailed investigation of subgradients, we refer the reader to one of the works of 
e.g. [Roc97, BV04, Nem05, KZ05, BL06]. 

2.2 A Duality for Convex Optimization over Compact Domain 

For a given point x & D, and any choice of a subgradient dx G df{x), we define a dual function 
value 

uj{x,dx) := Toain f{x) + {y - x,dx) . (3) 

In other words uj{x,dx) G M is the minimum of the linear approximation to / defined by the 
subgradient dx at the supporting point x, where the minimum is taken over the domain D. This 
minimum is always attained, since D is compact, and the linear function is continuous in y. 

By the definition of the subgradient — as lying below the graph of the function / — we readily 
attain the property of weak-duality, which is at the core of the optimization approach we will 
study below. 

Lemma 1 (Weak duality). For all pairs x,y £ D, it holds that 

^{x,dx) < f{y) 

Proof. Immediately from the definition of the dual a;(., .): 

Lj{x,dx)= min^gz) /(x) + (z - 

< f{x) + {y - x,dx) 

< f{y)- 

Here the last inequality is by the definition (2) of a subgradient. □ 

Geometrically, this fact can be understood as that any function value /(y), which is "part of" 
the graph of /, always lies higher than the minimum over any linear approximation (given by 
dx) to /. 

In the case that / is differentiable, there is only one possible choice for a subgradient, namely 
dx = Vf{x), and so we will then denote the (unique) dual value for each x by 

uj{x) := uj{x, V/(x)) = min f{x) + {y - x, V/(x)) . (4) 

■^The assumption that our objective function / is lower semicontinuous on D, is equivalent to the fact that its 
epigraph — i.e. the set {{x,t) G D x R | f > f{x)} of all points lying on or above the graph of the function 
— is closed, see also [KZ05, Theorem 7.1.2]. 
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The Duality Gap as a Measure of Approximation Quality. The above duality concept allows 
us to compute a very simple measure of approximation quality, for any candidate solution x £ D 
to problem (1). This measure will be easy to compute even if the true optimum value f{x*) is 
unknown, which will very often be the case in practical applications. The duality gap g{x,dx) 
at any point x £ D and any subgradient subgradient dx G df{x) is defined as 

g{x,dx) := f{x) - uj{x,dx) = max {x - y , d^) , (5) 

or in other words as the difference of the current function value f{x) to the minimum value 
of the corresponding linearization at x, taken over the domain D. The quantity g{x,dx) = 
f{x) — uj{x, dx) will be called the duality gap at x, for the chosen dx- 

By the weak duality Lemma 1, we obtain that for any optimal solution x* to problem (1), it 
holds that 

9{x, dx) > fix) - fix*) > yxeD, Wx G dfix) . (6) 

Here the quantity fix) — fix*) is what we call the primal error at point x, which is usually 
impossible to compute due to x* being unknown. The above inequality (6) now gives us that 
the duality gap — which is easy to compute, given dx — is always an upper bound on the 
primal error. This property makes the duality gap an extremely useful measure for example as 
a stopping criterion in practical optimizers or heuristics. 

We call a point x €z X an e- approximation if gix,dx) < e for some choice of subgradient 
dx G dfix). 

For the special case that / is differentiable, we will again use the simpler notation gix) for 
the (unique) duality gap for each x, i.e. 

gix) := gix, V/(x)) = max {x - y, V/(x)) . 

y£D 

Relation to Duality of Norms. In the special case when the optimization domain D is given 
by the unit ball of some norm on the space X, we observe the following: 

Observation 2. For optimization over any domain D = {x £ X \ ||x||<l} being the unit ball 
of some norm \\.\\, the duality gap for the optimization problem min/(x) is given by 

gix,dx) = \\dx\\^ + {x,dx) , 
where \\.\\^ is the dual norm o/||.||. 

Proof. Directly by the definitions of the dual norm ||x||^ = sup||j^||<i(y, x), and the duality gap 
gix,dx) = maXj,gD {y,-dx) + {x,dx) as in (5). □ 



3 A Projection-Free First-Order Method for Convex Optimization 
3.1 The Algorithm 

In the following we will generalize the sparse greedy algorithm of [FW56] and its analysis 
by [ClalO] to convex optimization over arbitrary compact convex sets D C A' of a vector space. 
More formally, we assume that the space X is a Hilbert space, and consider problems of the 
form (1), i.e., 

minimize fix) . 

Here we suppose that the objective function / is differentiable over the domain D, and that for 
any x £ D, we are given the gradient V/(x) via an oracle. 
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The existing algorithms so far did only apply to convex optimization over the simplex (or 
convex hulls in some cases) [ClalO], or over the spectahedron of PSD matrices [Haz08], or then 
did not provide guarantees on the duality gap. Inspired by the work of Hazan [Haz08], we 
can also relax the requirement of exactly solving the linearized problem in each step, to just 
computing approximations thereof, while keeping the same convergence results. This allows for 
more efficient steps in many applications. 

Also, our algorithm variant works for arbitrary starting points, without needing to compute 
the best initial "starting vertex" of D as in [ClalO]. 

The Primal-Dual Idea. We are motivated by the geometric interpretation of the "poor man's" 
duality gap, as explained in the previous Section 2. This duality gap is the maximum difference 
of the current value f{x), to the linear approximation to / at the currently fixed point x, where 
the linear maximum is taken over the domain D. This observation leads to the algorithmic idea 
of directly using the current linear approximation (over the convex domain D) to obtain a good 
direction for the next step, automatically staying in the feasible region. 

The general optimization method with its two precision variants is given in Algorithm 1. 
For the approximate variant, the constant C/ > is an upper bound on the curvature of the 
objective function /, which we will explain below in more details. 

Algorithm 1 Greedy on a Convex Set 

Input: Convex function /, convex set D, target accuracy e 
Output: e-approximate solution for problem (1) 
Pick an arbitrary starting point x^^^ £ D 
for = . . . oo do 
Let a := ^2 

Compute s := ExaCtLineaR (V/(a;(^')), D) {Solve the linearized primitive task exactly} 
— or — 

Compute s := ApproxLineaR (V/(a;(*'')), D, aCf) {Approximate the Hnearized problem} 
Update := x'^^^ + a{s - x^^^) 

end for 



The Linearized Primitive. The internal "step direction" procedure ExactLinear(c, D) used 
in Algorithm 1 is a method that minimizes the linear function (x, c) over the compact convex 
domain D. Formally it must return a point s € D such that (s, c) = min(y,c). In terms of 

the smooth convex optimization literature, the vectors y that have negative scalar product with 
the gradient, i.e. (y, V/(x)) < 0, are called descent directions, see e.g. [BV04, Chapter 9]. The 
main difference to classical convex optimization is that we always choose descent steps staying 
in the domain D, where traditional gradient descend techniques usually use arbitrary directions 
and need to project back onto D after each step. We will comment more on this analogy in 
Section 3.7. 

In the alternative interpretation of our duality concept from Section 2, the linearized sub-task 
means that we search for a point s that "realizes" the current duality gap g{x), that is the 
distance to the linear approximation, as given in (5). 

In the special case that the set D is given by an intersection of linear constraints, this sub-task 
is exactly equivalent to linear programming, as already observed by [FW56, Section 6]. However, 
for many other representations of important specific feasible domains D, this internal primitive 
operation is significantly easier to solve, as we will see in the later sections. 
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Figure 1: Visualization of a step of Algorithm 1, moving from the current point x = a;*-'^-* towards a linear 
niinimizer s £ D. Here the two-dimensional domain D is part of the ground plane, and we plot 
the function values vertically. Visualization by Robert Carnecky. 

The alternative approximate variant of Algorithm 1 uses the procedure ApproxLinear (c, D, e') 
as the internal "step-direction generator". Analogously to the exact primitive, this procedure 
approximates the minimum of the linear function {x, c) over the convex domain D, with ad- 
ditive error e' > 0. Formally ApproxLinear (c, D, e') must return a point s £ D such that 

(s, c) < min {y, c) + e'. For several applications, this can be done significantly more efficiently 
j/eD 

than the exact variant, see e.g. the applications for semidefinite programming in Section 7. 

The Curvature. Everything we need for the analysis of Algorithm 1 is that the linear approx- 
imation to our (convex) function / at any point x does not deviate from / by too much, when 
taken over the whole optimization domain D. 

The curvature constant Cj of a convex and differentiable function / : M" — )■ M, with respect 
to the compact domain D is defined as. 

Cf := sup ^ (/(y) - fix) - {y - x, Vf{x))) . (7) 

x,sGD, 

ae[o,i], 

y=x+a(s—x) 

A motivation to consider this quantity follows if we imagine our optimization procedure at the 
current point x = x^^\ and choosing the next iterate as y = x'-'^"''^^ := x + a{s — x). Bounded 
Cf then means that the deviation of / at y from the "best" linear prediction given by V/(x) is 
bounded, where the acceptable deviation is weighted by the inverse of the squared step-size a. 
For linear functions / for example, it holds that C/ = 0. 

The defining term f{y) — f{x) — {y — x,dx) is also widely known as the Bregman divergence 
defined by /. The quantity Cf turns out to be small for many relevant applications, some of 
which we will discuss later, or see also [ClalO]. 

The assumption of bounded curvature C/ is indeed very similar to a Lipschitz assumption on 
the gradient of /, see also the discussion in Sections 3.4 and 3.7. In the optimization literature, 
this property is sometimes also called Cf -strong smoothness. 

It will not always be possible to compute the constant Cf exactly. However, for all algorithms 
in the following, and also their analysis, it is sufficient to just use some upper bound on Cf. We 
will comment in some more details on the curvature measure in Section 3.4. 
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Convergence. The following theorem shows that after O(^) many iterations, Algorithm 1 
obtains an e-approximate solution. The analysis essentially follows the approach of [ClalO], 
inspired by the earlier work of [FW56, DH78, Pat93] and [Zlia03]. Later, in Section 4.2, we 
will show that this convergence rate is indeed best possible for this type of algorithm, when 
considering optimization over the unit-simplex. More precisely, we will show that the dependence 
of the sparsity on the approximation quality, as given by the algorithm here, is best possible up 
to a constant factor. Analogously, for the case of semidefinite optimization with bounded trace, 
we will prove in Section 7.4 that the obtained (low) rank of the approximations given by this 
algorithm is optimal, for the given approximation quality. 

Theorem 3 (Primal Convergence). For each k>l, the iterate x^''^ of the exact variant of Algo- 
rithm 1 satisfies 

/(x('=))-/(x*)<^, 

where x* G D is an optimal solution to problem (1). For the approximate variant of Algorithm 1, 
it holds that 

/(^W)_/(^*)<i^ . 

(In other words both algorithm variants deliver a solution of primal error at most e after O(^) 
many iterations.) 

The proof of the above theorem on the convergence-rate of the primal error crucially depends 
on the following Lemma 4 on the improvement in each iteration. We recall from Section 2 
that the duality gap for the general convex problem (1) over the domain D is given by g{x) = 
max {x — s, V/(x)). 

Lemma 4. For any step x^'''^^^ := x^^'^ + a{s — x^^^) with arbitrary step-size a E [0, 1], it holds 
that 

/(x^'^+i)) < /(x^'^)) - ag{x^^^) + a^Cj 

if s is given by s := ExactLinear (V/(x), D). 

// the approximate primitive s := AppROxLiNEAR (V/(x), D, aC/) is used instead, then it 
holds that 

/(xC^+i)) < /(xW) - a5(x('=)) + 2a^Cf . 

Proof. We write x := x^^\ y := x^^^^^ = x-|-a(s — x), and dx ■= V/(x) to simplify the notation, 
and first prove the second part of the lemma. We use the definition of the curvature constant 
Cfol our convex function /, to obtain 

f{y)= f{x + a{s-x)) 

< f{x)-^a{s — x,dx)-^a'^Cf. 

Now we use that the choice of s := ApproxLinear (d^^, D, e') is a good "descent direction" on 

the linear approximation to / at x. Formally, we are given a point s that satisfies {s, dx) < 

min(y, dx) + e', or in other words we have 
yen 

{s - X, dx) < minygD(y, dx) - (x, dx) + e' 
= -g{x,dx) + e' , 

from the definition (5) of the duality gap g{x) = g{x,dx)- Altogether, we obtain 

f{y)< fix) + a{-g{x)+e') + a^Cf 
= f{x)-ag{x) + 2a^Cf , 
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the last equality following by our choice of e' = ctCf. This proves the lemma for the approximate 
case. The first claim for the exact linear primitive ExactLinear() follows by the same proof 
for e' = 0. □ 

Having Lemma 4 at hand, the proof of our above primal convergence Theorem 3 now follows 
along the same idea as in [ClalO, Theorem 2.3] or [Haz08, Theorem 1]. Note that a weaker 
variant of Lemma 4 was already proven by [FW56]. 

Proof of Theorem 3. From Lemma 4 we know that for every step of the exact variant of Algo- 
rithm 1, it holds that < /(x^^)) - ac/(x(^)) + a'^Cf. 

Writing h{x) := f{x) — f{x*) for the (unknown) primal error at any point x, this implies that 

< - ag{x^''^) + a'^Cf 

< h{x^^'^)-ah{x^^'^) + a'^Cf (8) 
= {l-a)h{x^^^) + a^Cf , 

where we have used weak duality h{x) < g{x) as given by in (6). We will now use induction 
over k in order to prove our claimed bound, i.e., 

h{x^k+i)^ < A: = 0...oo. 



The base-case A: = follows from (8) applied for the first step of the algorithm, using a = a^"^ = 
^ = 1 

0+2 

Now considering /c > 1, the bound (8) gives us 

h{x^k+i)) < (1 _ aW)/i(xW) + a('=)^C/ 

< (l-42)a + (42m> 

where in the last inequality we have used the induction hypothesis for x^^\ Simply rearranging 
the terms gives 



4C,' ^ 1 



\k+2 (fc+2)2 
^ 4C/ fc+2-1 
fc+2 fc+2 
4CV fc+2 
— fe+2 fc+3 

which is our claimed bound for A; > 1. 

The analogous claim for Algorithm 1 using the approximate linear primitive ApproxLinear() 
follows from the exactly same argument, by replacing every occurrence of in the proof here 
by 2Cf instead (compare to Lemma 4 also). □ 

3.2 Obtaining a Guaranteed Small Duality Gap 

From the above Theorem 3 on the convergence of Algorithm 1, we have obtained small primal 
error. However, the optimum value f{x*) is unknown in most practical applications, where we 
would prefer to have an easily computable measure of the current approximation quality, for 
example as a stopping criterion of our optimizer in the case that C/ is unknown. The duality 
gap g{x) that we defined in Section 2 satisfies these requirements, and always upper bounds the 
primal error f{x) — f{x*). 
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By a nice argument of Clarkson [ClalO, Theorem 2.3], the convergence on the simplex opti- 
mization domain can be extended to obtain the stronger property of guaranteed small duality 
gap ^(x^'^)) < e, after at most O(^) many iterations. This stronger convergence result was not 
yet known in earlier papers of [FW56, DH78, Jon92, Pat93] and [ZhaOS]. Here we will generalize 
the primal-dual convergence to arbitrary compact convex domains. The proof of our theorem 
below again relies on Lemma 4. 



Theorem 5 (Primal-Dual Convergence). Let K := . We run the exact variant of Algo- 

rithm 1 for K iterations (recall that the step-sizes are given by a^''^ := -j^, < k < K), 
and then continue for another K + 1 iterations, now with the fixed step-size a^^^ := -j^j^ for 
K < k <2K + 1. 

Then the algorithm has an iterate x^^\ K <k < 2K + 1, with duality gap hounded by 

< e . 

The same statement holds for the approximate variant of Algorithm 1, when setting K := 



instead. 



Proof. The proof follows the idea of [ClalO, Section 7]. 

By our previous Theorem 3 we already know that the primal error satisfies 
f{x*) < after K iterations. In the subsequent K -\- 1 iterations, we will now suppose that 

g{x^''^) always stays larger than We will try to derive a contradiction to this assumption. 

Putting the assumption g{x^^^) > into the step improvement bound given by Lemma 4, 
we get that 

< -aWia + ae^)^^, 

holds for any step size a^'^^ £ (0, 1]. Now using the fixed step-size a^'^^ = in the iterations 
k > K oi Algorithm 1, this reads as 

J{-^ ) jy-^ )^ K+2 K+2 ^ (K+2)2 '^f 

— 4Cj- 

- W+W 

Summing up over the additional steps, we obtain 

2K+1 

^ {-^ ^ ^) (ii-+2)2 — K+2 ' 

which together with our known primal approximation error -/(x*) < ^ would result 

in /(x^^^"*"^^) — f{x*) < 0, a contradiction. Therefore there must exist k, K < k < 2K -\- 1, with 

The analysis for the approximate variant of Algorithm 1 follows using the analogous second 
bound from Lemma 4. □ 

Following [ClalO, Theorem 2.3], one can also prove a similar primal-dual convergence theorem 
for the line-search algorithm variant that uses the optimal step-size in each iteration, as we will 
discuss in the next Section 3.3. This is somewhat expected as the line-search algorithm in each 
step is at least as good as the fixed step-size variant we consider here. 
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3.3 Choosing the Optimal Step-Size by Line-Search 

Alternatively, instead of the fixed step-size a = -^^^ ™- Algorithm 1 , one can also find the optimal 
a G [0, 1] by line-search. This will not improve the theoretical convergence guarantee, but might 
still be considered in practical applications if the best a is easy to compute. Experimentally, we 
observed that line-search can improve the numerical stability in particular if approximate step 
directions are used, which we will discuss e.g. for semidefinite matrix completion problems in 
Section 11.5. 

Formally, given the chosen direction s, we then search for the a of best improvement in the 
objective function /, that is 

Q := argmin / (x^'') + a(s - x^''))) . (9) 
ae[o,i] ^ ^ 

The resulting modified version of Algorithm 1 is depicted again in Algorithm 2, and was precisely 
analyzed in [ClalO] for the case of optimizing over the simplex. 

Algorithm 2 Greedy on a Convex Set, using Line-Search 
Input: Convex function /, convex set D, target accuracy e 
Output: e-approximate solution for problem (16) 
Pick an arbitrary starting point a;^*^^ G D 
for A; = . . . oo do 

Compute s := ExactLinear (yf{x^^^),D) 
— or — 

Compute s := ApproxLinear (v /{x'^'^^), D, 

Find the optimal step-size a := argmin / (x^''^ + a(s — x^''^)] 

ae[o,i] ^ ^ 

Update x^'^+i) := x'-''^ + a{s - x^^^) 

end for 



In many cases, we can solve this line-search analytically in a straightforward manner, by differ- 
entiating the above expression with respect to a: Consider fa ■= f (^(a)"^^) — f (^^^^ + a (s — x^^^)) 
and compute 




If this equation can be solved for a, then the optimal such a can directly be used as the step-size 
in Algorithm 1, and the convergence guarantee of Theorem 3 still holds. This is because the 
improvement in each step will be at least as large as if we were using the older (potentially 
sub-optimal) fixed choice of a = -j^. Here we have assumed that a^''^ G [0,1] always holds, 
which can be done when using some caution, see also [ClalO]. 

Note that the line-search can also be used for the approximate variant of Algorithm 1, where 
we keep the accuracy for the internal primitive method ApproxLinear (\7 f{x^^^), D,e') fixed 
to e' = aflxcdC*/ = Theorem 3 then holds as as in the original case. 

3.4 The Curvature Measure of a Convex Function 

For the case of differentiable / over the space X = M", we recall the definition of the curvature 
constant Cf with respect to the domain D C M", as stated in (7), 

Cf := sup ^ (/(y) - fix) - (y - xfVf{x)) . 

Q6[0,l], 
y=x+a{s—x) 
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An overview of values of for several classes of functions / over domains that are related to 
the unit simplex can be found in [ClalO] . 

Asymptotic Curvature. As Algorithm 1 converges towards some optimal solution x* , it also 
makes sense to consider the asymptotic curvature Cj, defined as 

C}:= sup J^(/(y)-/(x*)-(y-x*)^V/(x*)) . (11) 
ae[o,i], 

y=x* +a{s—x*) 

Clearly Cj < Cj. As described in [ClalO, Section 4.4], we expect that as the algorithm converges 
towards x* , also the improvement bound as given by Lemma 4 should be determined by Cj+o{l) 
instead of C/, resulting in a better convergence speed constant than given Theorem 3, at least 
for large k. The class of strongly convex functions is an example for which the convergence of 
the relevant constant towards is easy to see, since for these functions, convergence in the 
function value also imlies convergence in the domain, towards a unique point x* , see e.g. [BV04, 
Section 9.1.2]. 

Relating the Curvature to the Hessian Matrix. Before we can compare the assumption of 
bounded curvature C/ to a Lipschitz assumption on the gradient of /, we will need to relate Cf 
to the Hessian matrix (matrix of all second derivatives) of /. 

Here the idea described in [ClalO, Section 4.1] is to make use of the degree-2 Taylor-expansion 
of our function / at the fixed point function of a, which is 

2 

/(^ + a{s - x)) = fix) + a{s - xfVf{x) + ^{s - xfV^f{z){s - x) , 

where z is a point on the line-segment [x, y] C D C M.'^ between the two points x G and 
y = X + a{s — x) £ M". To upper bound the curvature measure, we can now directly plug in 
this expression for f{y) into the above definition of Cj, obtaining 

Cf< sup hy-xfV^f{z){y-x) . (12) 

x,'y£D, ^ 
zGlx,y]CD 

From this bound, it follows that C/ is upper bounded by the largest eigenvalue of the Hessian 
matrix of /, scaled with the domain's Euclidean diameter, or formally 

Lemma 6. For any twice differentiable convex function f over a compact convex domain D, it 
holds that 

Cf <\ diam(Z))2 . sup A,nax (V2/(^)) • 

Note that as / is convex, the Hessian matrix V^/(z) is positive semidefinite for all z, see 
e.g. [KZ05, Theorem 7.3.6]. 

Proof. Applying the Cauchy-Schwarz inequality to (12) for any x,y £ D (as in the definition of 
Cf), we get 

{y - x)'^V'^f{z){y -x) < \\y - x\\^ \\V^f{z){y - x)\\^ 

< ||y-.||^||vV(z)||^^^^ 

< diam(i?)2 . sup A^ax (v2/(^)) • 

The middle inequality follows from the variational characterization of the matrix spectral norm, 

\\Ax\\ 

i.e. Il^ll^pgc := sup^j-^o \\x\\ • Finally, in the last inequality we have used that by convexity of 
/, the Hessian matrix V^/(2;) is PSD, so that its spectral norm is its largest eigenvalue. □ 
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Note that in the case of D being the unit simplex (see also the following Section 4), we have 

that diam(A„) = \/2, meaning the scaling factor disappears, i.e. C/ < sup Amax (^^/(z)). 

zeA„ 

Bounded Curvature vs. Lipschitz-Continuous Gradient. Our core assumption on the given 
optimization problems is that that the curvature C/ of the function is bounded over the domain. 
Equivalently, this means that the function does not deviate from its linear approximation by too 
much. Here we will explain that this assumption is in fact very close to the natural assumption 
that the gradient V/ is Lipschitz-continuous, which is often assumed in classical convex opti- 
mization literature, where it is sometimes called Cf -strong smoothness, see e.g. [Ncm05, KSST09] 
(or [d'A08] if the gradient information is only approximate). 

Lemma 7. Let f be a convex and twice differentiable function, and assume that the gradient V/ 
is Lipschitz-continuous over the domain D with Lip schitz- constant L > 0. Then 

Cf<^ diam(L>)2L . 

Proof. Having ||V/(y) — V/(x)||2 < L||y — x||2 Vx,y € -D by the Cauchy-Schwarz inequality 
implies that (y — x)'^{Vf{y) — V/(x)) < L\\y — xHg, so that 

f{y) < fix) + {y- xfVfix) + ^\\y- x\\l . (13) 

If / is twice differentiable, it can directly be seen that the above condition implies that 
L-Ih V2/(z) holds for the Hessian of /, that is Amax (V2/(-z)) < L. 

Together with our result from Lemma 6, the claim follows. □ 

The above bound (13) which is implied by Lipschitz-continuous gradient means that the 
function is not "curved" by more than L in some sense, which is an interesting property. In fact 
this is exactly the opposite inequality compared to the property of strong convexity, which is 
an assumption on the function / that we do not impose here. Strong convexity on a compact 
domain means that the function is always curved at least by some constant (as our L). We just 
note that for strongly convex functions, "accelerated" algorithms with an even faster convergence 
of (meaning O(^) steps) do exist [Nes04, Nes07a]. 

3.5 Optimizing over Convex Hulls 

In the case that the optimization domain D is given as the convex hull of a (finite or infinite) 
subset V d X , i.e. 

D = conv(y) , 

then it is particularly easy to solve the linear optimization subproblems as needed in our Algo- 
rithm 1. Recall that conv(y) is defined as the set of all finite convex combinations Oi^'i for 
a finite subset {vi, . . . , Vk} C V , while V can also be an infinite set. 

Lemma 8 (Linear Optimization over Convex Hulls). Let D = conv(y) for any subset V <Z X , and 
D compact. Then any linear function y i— t- {y,c) will attain its minimum and maximum over D 
at some "vertex" v & V . 

Proof. W.l.g. we will only show the case for the maximum. Let s € -D be a point attaining the 
linear optimum (s,c) = maxyg£)(y, c). Then by definition of D, we have that s = ajVj, 
meaning that s is the weighted average of some finite set of "vertices" vi,...,Vk G V, with 
oti > 0, Ylii Oil = 1. By linearity of the inner product. 
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(s, c) = ( ^ aiVi, c) =^ ai{vi,c) , 

\i=l I i=l 

and therefore we must have that {vi,c) > (s,c) for at least one of the indices i, meaning that 
Vi GV is also attaining the linear maximum. □ 

In the following we will discuss several application where this simple fact will be useful to 
solve the linearized subproblems ExactLinear() more efficiently, as the set V is often much 
easier to describe than the full compact domain D. 

The setting of convex optimization over a convex hull in a vector space was already studied 
by [Zha03]. There, each iteration of the optimizer consists of greedily selecting the point (or 
"vertex" ) of V which promises best improvement. [Zha03] then gave a similar primal convergence 
guarantee as in our Theorem 3 (but no primal-dual convergence result on general convex hulls was 
known so far, to the best of our knowledge). The above Lemma 8 in a sense explains the relation 
to our linearized internal problem. The main difference is that the algorithm of [Zha03] always 
evaluates the original non- linear function / at all vertices V, while our slightly more general 
framework only relies on the linear subproblem, and allows for arbitrary means to approximately 
solve the subproblem. 

3.6 Randomized Variants, and Stochastic Optimization 

For a variety of classes of our convex optimization problems, randomization can help to solve 
the linearized subproblem more efficiently. This idea is strongly related to online and stochastic 
optimization, see e.g. [Ncsll], and also the popular stochastic gradient descent (SGD) tech- 
niques [BotlO]. 

We can also apply such cheaper randomized steps in our described framework, instead of 
deterministically solving the internal linear problem in each iteration. Assumed that the user of 
our method is able to decompose the linearized problem in some arbitrary way using random- 
ization, and if the randomization is such that the linearized problem will be solved "accurately 
enough" with some probability p > in each iteration, then our convergence analysis still holds 
also in this probabilistic setting as follows: 

Formally, we assume that we are given access to a randomized procedure RandomLinear (c, D, e'), 
which returns a point s £ D such that (s, c) < min (y, c) + e' with probability at least p > 0. In 

j/GD 

other words, with a positive probabihty, RandomLinear() wih behave like ApproxLinear(). 
In each iteration of the line-search variant of our algorithm (see Algorithm 2), we will now use 
that randomized internal procedure instead. The expected improvement given by a step towards 
s = RandomLinear() is at least p times the amount given in Lemma 4. (Here we have used 
that in the events of "failure" of RandomLinear(), the objective function value will at least 
not become worse, due to the use of line-search). 

In other words if we perform ^ times more iterations than required for the deterministic 
Algorithm 1, then we have that the convergence by Theorem 3 also holds for the randomized 
variant described here. 

Stochastic Gradient Descent (SGD). A classical example is when the linearized problem is 
given by simply finding the maximum over say n coordinates, as we will e.g. see in the following 
Sections 4 and 5 for optimizing over the simplex, or over bounded £i-norm. In this case, by 
sampling uniformly at random, with probability ^ we will pick the correct coordinate, for which 
the step improvement is as in the deterministic Algorithm L Therefore we have obtained the 
same convergence guarantee as for the deterministic algorithm, but the necessary number of 
steps is multiplied by n. 
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For unconstrained convex optimization, the convergence of SGD and other related methods 
was analyzed e.g. in [Ncsll] and also the earlier paper [NcsOTb, Section 6]. Also here, a 
comparable slow-down was observed when using the cheaper randomized steps. 



3.7 Relation to Classical Convex Optimization 

Relation to Gradient Descent and Steepest Descent. The internal linear optimizer in our 
Algorithm 1 can also be interpreted in terms of descent directions. Recall that all vectors y 
that have negative scalar product with the current gradient, i.e. {y,Vf{x)) < 0, are called 
descent directions, see e.g. [BV04, Chapter 9.4]. Also observe that (y, V/(x)) is the directional 
derivative of / in direction of y if y is of unit length. Our method therefore chooses the best 
descent direction over the entire domain D, where the quality is measured as the best possible 
absolute improvement as suggested by the linearization at the current point x. In any iteration, 
this will crucially depend on the global shape of the domain D, even if the actual step-size a^^^ 
might be very small. 

This crucially contrasts classical gradient descent techniques, which only use local information 
to determine the step-directions, facing the risk of walking out of the domain D and therefore 
requiring projection steps after each iteration. 



Relation to Inaccurate and Missing Gradient Information. The ability of our Algorithm 1 
to deal with only approximate internal linear optimizers as in ApproxLinear() is also related 
to existing methods that assume that gradient information is only available with noise, or in a 
stochastic or sampling setting. 

For the case of optimizing smooth convex functions, [d'A08] has used a similar measure of 
error, namely that the linearization given by the "noisy" version dx of the gradient V/(a;) does 
not differ by more than say e' when measured over the entire domain D, or formally 



{y - z,dx) - {y- z,Vf{x)) 



<e', yx,y,z e D . 



This assumption is similar, but stronger than the additive approximation quality that we require 
in our above setting (we only need that the linearized optimal values are within e'). Also, the 
algorithm as well as the analysis in [d'A08] are more complicated than the method proposed 
here, due to the need of projections and proximal operators. 

We have discussed the case where gradient information is available only in a stochastic oracle 
(e.g. such that the gradient is obtained in expectation) in the above Subsection 3.6. For an 
overview of related randomized methods in unconstrained convex optimization, we refer the 
reader to the recent work by [Ncsll], which also applies when the gradient itself is not available 
and has to be estimated by oracle calls to the function alone. 

If gradient information can be constructed in any way such that the linearized problem 
ApproxLinear() can be solved to the desired additive error, then our above analysis of Algo- 
rithm 1 will still hold. 



Relation to Mirror Descent, Proximal Methods and Conjugate Functions. Our proposed 
method is related to mirror descent as well as proximal methods in convex optimization, but our 
approach is usually simpler. The mirror descent technique originates from e.g. [BT03, BTN05]. 
For a brief overview of proximal methods with applications to some of the classes of sparse 
convex optimization problems as studied here, we refer to [BJMOll, Section 3]. 

To investigate the connection, we write fun\xiy) ■= f{x) + {y — x,dx) for the linearization 
given by the (sub)gradient dx = V/(x) at a fixed point x G D. A variant of mirror descent. 
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see e.g. [BTN05, Hazll] is to find the next iterate y as the point maximizing the Bregman 
divergence 

fiy) - fix) -{y-x, 4) = /(y) - A„|,(2/) (14) 

relative to the currently fixed old iterate x. This is the same task as maximizing the gap between 
the function f(y) and its linear approximation at x, or equivalently we evaluate the conjugate 
function f*{z) := sup(y, z) — f{y) at z = dx- The definition of the conjugate dual is also known 

y€D 

as Fenchel duality, see e.g. [BL06]. In [Nem05], the conjugate function is also called the Legendre 
transformation. 

However in our approach, the inner task ExACTLiNEAR((i2:, D) as well as AppROxLiNEAR((i2:, D, e') 
is a simpler linear problem. Namely, we directly minimize the linearization at the current point 
X, i.e. we maximize 

- fix) - {y-x, dx) = -fnn\x(.y) (15) 

and then move towards an approximate maximizer y. In terms of Fenchel duality, this simpler 
linear problem is the evaluation of the conjugate dual of the characteristic function of our domain 
D, i.e. 

l}){z) := sui>{y,z) - 1d{v) , 
ye AT 

where this function is evaluated at the current subgradient z = dx- The characteristic function 
\d : — 7- M of a set D C is defined as lz)(y) = for y & D and loiy) = oo otherwise. 

Compared to our algorithm, mirror descent schemes require a "projection" step in each it- 
eration, sometimes also called the proximal or mirror operator. This refers to minimizing the 
linearization plus a strongly convex prox-function that punishes the distance to the starting 
point. If the squared Euclidean norm is used, the mirror operator corresponds to the standard 
projection back onto the domain D. Our method uses no such prox-function, and neither is 
the zero-function a strongly convex one, as would be required for mirror descent to work. It is 
expected that the computational cost per iteration of our method will in most application cases 
be lower compared to mirror descent schemes. 

For convex optimization over the simplex, which we will study in more details in the following 
Section 4, [BT03] have proposed a mirror descent algorithm, obtaining a convergence of f{x^^'^) — 
f{x*) < \^2lnn^. This however is worse than the convergence of our methods as given by 

Theorem 3. Our convergence is independent of the dimension n, and goes with ^ instead of 
Also the earlier paper by [BTMNOl] only obtained a convergence of O(^) for the case of 
Lipschitz-continuous convex functions over the simplex. 

The NERML optimizer by [BTN05] is a variant of mirror descent that memorizes several past 
linearizations of the objective function, to improve the available knowledge about the current 
function landscape. It is an open research question if this idea could also help in our setting 
here, or for stochastic gradient descent schemes [BotlO]. 

4 Sparse Approximation over the Simplex 

As a first application of the above general scheme, we will now consider optimization problems 
defined over the unit simplex, or in other words the non- negative vectors of ^i-norm equal to 
one. This will serve as a warm-up case before considering ^i-norm regularized problems in the 
next Section 5. 

Our approach here will allow us to understand the best achievable sparsity of approximate 
solutions, as a function of the approximation quality, as already shown by [ClalO]. 

In particular, we will show that our main Algorithm 1 on page 8 and its analysis do lead to 
Clarkson's approach [ClalO] for optimizing over the simplex. In this case, it was already known 



18 



that sparsity O(^) can always be achieved by applying Algorithm 1 to the simplex domain, 
see [ClalO]. We will also show that this is indeed optimal, by providing an asymptotically 
matching lower bound in Section 4.2. Also, our analysis holds even if the linear subproblems 
are only solved approximately, and allows arbitrary starting points, in contrast to [ClalO]. 

Having this efficient algorithm giving sparsity O(^) is in particularly attractive in view of the 
computational complexity of vector cardinality minimization, which is known to be NP-hard, by 
a reduction to Exact-Cover, see [Nat95]. Vector cardinality minimization here refers to finding 
the sparsest vector that is an e-approximation to some given convex minimization problem. 
Formally, finding the sparsest x that satisfies \\Ax — ftjlg < e for given A, b and e. 

Set-Up. We suppose that a basis has been chosen in the space X, so that we can assume 
X = with the standard inner product {x,y) = x'^y. Here we consider one special class of the 
general optimization problems (1), namely we optimize over non- negative vectors that sum up 
to one, that is 

minimize fix) 

s.t. \\x\\^ = l, (16) 
X > . 

In the following we write A„ := {x € M" |x>0, ||x||-|^ = l} for the unit simplex in R". As 
the objective function / is now defined over M", all subgradients or gradients of / will also be 
represented by vectors in M" in the following. 

Note that the alternative case of optimizing under an inequality constraint ||x||^ < 1 instead 
of = 1 can easily be reduced to the above form (16) by introducing a new "slack" variable. 
Formally, one uses vectors x = (xi, . . . , x„, Xn+i) G M"+-'^ instead and optimizes the function 
/(x) := /(xi, . . . ,x„) over the simplex domain ||x||-|^ = 1, x > 0. 

Coresets. The coreset concept was originally introduced in computational geometry by [BHPI02] 
and [APV02]. For point-set problems, the coreset idea refers to identifying a very small subset 
(coreset) of the points, such that the solution just on the coreset is guaranteed to be a good 
approximation to original problem. Here for general convex optimization, the role of the coreset 
points is taken by the non-zero coordinates of our sought vector x instead. The coreset size then 
corresponds to the sparsity of x. 

Formally if there exists an e-approximate solution x €z D CI M" to the convex optimiza- 
tion problem (1), using only k many non-zero coordinates, then we say that the corresponding 
coordinates of x form an e-coreset of size k for problem (1). 

In other words, the following upper and lower bounds of 0(1) on the sparsity of approxi- 
mations for problem (16) are indeed matching upper and lower bounds on the coreset size for 
convex optimization over the simplex, analogous to what we have found in the geometric problem 
setting in [GJ09]. 

4.1 Upper Bound: Sparse Greedy on the Simplex 

Here we will show how the general algorithm and its analysis from the previous Section 3 do 
in particular lead to Clarkson's approach [ClalO] for minimizing any convex function over the 
unit simplex. The algorithm follows directly from Algorithm 1, and will have a running time 
of O (1) many gradient evaluations. We will crucially make use of the fact that every linear 
function attains its minimum at a vertex of the simplex A„. Formally, for any vector c G M", 
it holds that min s'^c = mincj . This property is easy to verify in the special case here, but is 

also a direct consequence of the small Lemma 8 which we have proven for general convex hulls. 
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if we accept that the unit simplex is the convex hull of the unit basis vectors. We have obtained 
that the internal linearized primitive can be solved exactly by choosing 

ExactLinear (c, A„) := with i = argmin Cj . 



Algorithm 3 Sparse Greedy on the Simplex 
Input: Convex function /, target accuracy e 
Output: e-approximate solution for problem (16) 
Set := ei 
for A; = . . . oo do 

Compute i := argmin^ {'V f{x^^^)^ . 

Let a := ^ 

Update := x^'^) + a{ei - x^^^) 

end for 

Observe that in each iteration, this algorithm only introduces at most one new non-zero 
coordinate, so that the sparsity of x^^^ is always upper bounded by the number of steps k, plus 
one, given that we start at a vertex. Since Algorithm 3 only moves in coordinate directions, it 
can be seen as a variant of coordinate descent. The convergence result directly follows from the 
general analysis we gave in the previous Section 3. 

Theorem 9 ([ClalO, Theorem 2.3], Convergence of Sparse Greedy on the Simplex). For each k > I, 
the iterate x^^'^ of Algorithm 3 satisfies 

j\ ) J\ J - ^_^2 
where x* G A„ is an optimal solution to problem (16). 



4C/ 

£ 



+ 1 = O (-) many steps , it has an iterate 



Furthermore, for any e > 0, after at most 2 
x^^) of sparsity O (^), satisfying g{x^^'^) < e. 

Proof. This is a corollary of Theorem 3 and Theorem 5. □ 

Duality Gap. We recall from Section 2 that the duality gap (5) at any point x £ A„ is easily 
computable from any subgradient, and in our case becomes 

g{x, dx) = x^dx - min {dx)i , and 

i 

g{x) = x'^V f{x) — min(V/(a;)) 



Here we have again used the observation that linear functions attain their minimum at a vertex 

of the domain, i.e, min c = mincj. 

seA„ i 

Applications. Many practically relevant optimization problems do fit into our setting (16) 
here, allowing the application of Algorithm 3. This includes linear classifiers such as support 
vector machines (SVMs), see also [GJ09], as well as kernel learning (finding the best convex 
combination among a set of base kernels) [BLJ04]. Some other applications that directly fit 
into our framework are ^2-support vector regression (SVR), AdaBoost [Zha03], mean- variance 
analysis in portfolio selection [Mar52], the smallest enclosing ball problem [BC07], and estimating 
mixtures of probability densities [ClalO]. For more applications we refer to [ClalO]. 



'Note that in order for our Theorem 5 on the bounded duahty gap to apply, the step-size in the second half 
of the iterations needs to be fixed to q''^' ~ ^qi^i see Section 3.2. This remark also applies to the later 
applications of our general Algorithm 1 that we discuss in the following. We already mentioned above that if 
line-search is used instead, then no such technicality is necessary, see also [ClalO]. 
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Line-Search for the Best Step-Size. In most applications it will be a straight-forward task to 
find the optimal step-size a G [0, 1] in each step instead, as described in Section 3.3. 

For the special case of polytope distance and SVM problems, the resulting method then exactly 
corresponds to Gilbert's geometric algorithm [Gil66], as shown in [GJ09]. Here the wording of 
"line-search" makes geometric sense in that we need to find the point s on a given line, such 
that s is closest to the origin. 



Away Steps. By performing more work with the currently non-zero coordinates, one can get 
the sparsity even smaller. More precisely the number of non-zeros can be improved close to 
instead of 2 as given by the above Theorem 9. The idea of away-steps introduced 

by [TY07] is to keep the total number of non-zero coordinates (i.e. the coreset size) fixed over 
all iterations, by removing the smallest non-zero coordinate from x after each adding step. For 
more background we refer to [ClalO, Algorithm 9.1]. 

4.2 Lower Bound on the Sparsity 

We will now show that sparsity O(^), as obtained by the greedy algorithm we analyzed in the 
previous section is indeed best possible, by providing a lower bound of In the language of 

coresets, this means we will provide a matching lower bound on the size of coresets for convex 
optimization over the simplex. Together with the upper bound, this therefore completely char- 
acterizes the trade-off between sparsity and approximation quality for the family of optimization 
problems of the form (16). The same matching sparsity upper and lower bounds will also hold 
for optimizing over the £i-ball instead, see Section 5. 

For the following lower bound construction we consider the differentiable function f{x) := 

1 1 1 1 2 — * 

This function has gradient V/(x) = 2x. Its curvature constant is Cf = 2, which 
follows directly from the definition (7), and the fact that here /(y) — f{x) — {y — x)^V/(x) = 
y'^y — x^x — {y — x)'^2x = \\x — y\\2, so that Cf = sup^, \\x — s\\2 = diam(A„)^ = 2. 

The following lemmata show that the sparse greedy algorithm of [ClalO] from Section 4.1 is 
indeed optimal for the approximation quality (primal as well as dual error respectively), giving 
best possible sparsity, up to a small multiplicative constant. 

Lemma 10. For f{x) := ||x||2, and 1 < k < n, it holds that 

min f(x) = — . 
xeA„ k 

card(x)<fc 

Proof. We prove the inequality min/(x) > i; by induction on k. 

X.. 

? any unit length vector x £ A„ having just a single non-zero entry, f{x) = 



Case 


k 


= 1 


2 ~ 


X 


1 ~ 


Case k > 1 



and write x = {1 — a)v + aei as the sum of two orthogonal vectors v and a unit basis vector ej, 
where v G A„ of sparsity < k — 1, Vi = 0, and a = Xi. So for every x £ A„ of sparsity < k, we 
therefore get that 

f{x) = \\x\\2 = x'^x 

= {{1 - a)v + aei)'^ {{1 - a)v + aei) 
= (1 - afv'^v + 

> (l-af^ + a^ 

> mino</3<i(l-/3)2^+/32 
_ 1 

~ fe- 
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In the first inequality we have applied the induction hypothesis for v G A„ of sparsity <k — \. 

Equality: The function value f{x) = ^ is obtained by setting k of the coordinates of x to ^ 
each. □ 

In other words for any vector x of sparsity card(x) = k, the primal error f{x) — f{x*) is always 
lower bounded by ^ — ^. For the duality gap g{x), the lower bound is even slightly higher: 

Lemma 11. For f{x) := ||x||2, and any k £ N, k < n, it holds that 

2 

g{x) > — Vx S A„ s.t. card(x) < k. 

Proof. g{x) = x^Vf{x) — minj(V/(x))j = 2{x'^x — minjXj). We now use miuj = because 
card(x) < n, and that by Lemma 10 we have x^x = f{x) > |. □ 

Note: We could also consider the function f(x) := 7 ||x||2 instead, for some 7 > 0. This / has 
curvature constant Cf = 27, and for this scaling, our above lower bound on the duality gap will 
also scale linearly, giving 

5 Sparse Approximation with Bounded £i-Norm 

In this second application case, will apply the general greedy approach from Section 3 in order 
to understand the best achievable sparsity for convex optimization under bounded £i-norm, as 
a function of the approximation quality. Here the situation is indeed extremely similar to the 
above Section 4 of optimizing over the simplex, and the resulting algorithm will again have a 
running time of O (^) many gradient evaluations. 

It is known that the vector ||.||j^-norm is the best convex approximation to the sparsity (car- 
dinality) of a vector, that is card(.). More precisely, the function ||.||-|^ is the convex envelope 
of the sparsity, meaning that it is the "largest" convex function that is upper bounded by the 
sparsity on the convex domain of vectors {x \ \\x\\^ < 1}. This can be seen by observing that 
card(x) > jj^jh , see e.g. [RFPIO]. We will discuss the analogous generalization to matrices in 

II II CXD 

the second part of this article, see Section 11, namely using the matrix nuclear norm as the 
"best" convex approximation of the matrix rank. 

Set- Up. Here we consider one special class of the general optimization problem (1), namely 
problems over vectors in M."' with bounded ||.||^-norm, that is 

minimize f(x) 

a;eK" (18) 
S.t. \\x\\^ < 1 . 

We write 0„ := {x G | ||x||-|^ < 1} for the ^i-ball in M". Note that one can simply rescale 
the function argument to allow for more general constraints ||x||^ < t for t > 0. Again we have 
X = M" with the standard inner product (x, y) = x'^y, so that also the subgradients or gradients 
of / are represented as vectors in M". 

The Linearized Problem. As already in the simplex case, the subproblem of optimizing a linear 
function over the ^i-ball is particularly easy to solve, allowing us to provide a fast implementation 
of the internal primitive procedure ExactLinear (c, {>„). 

Namely, it is again easy to see that every linear function attains its minimum/ maximum at a 
vertex of the ball ()n, as we have already seen for general convex hulls in our earlier Lemma 8, 
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and here {>„ = conv({ibej | i G [n]}). Here this crucial observation can also alternatively be 
interpreted as the known fact that the dual norm to the ^i-norm is in fact the £oo-norm, see also 
our earlier Observation 2. 

Observation 12. For any vector c G M", it holds that 

ej • sign(ci) E argmax y'^c 
y&On 

where i £ [n] is an index of a maximal coordinate of c measured in absolute value, or formally 
i £ arg maXj \cj\. 

Using this observation for c = — V/(x) in our general Algorithm 1, we therefore directly 
obtain the following simple method for £i-regularized convex optimization, as depicted in the 
Algorithm 4. 

Algorithm 4 Sparse Greedy on the £i-Ball 
Input: Convex function /, target accuracy e 
Output: e-approximate solution for problem (18) 
Set x(0) := 
for A; = . . . oo do 

Compute i := argmaXj |(V/(3;'^'^))) .| , 

and let s := Oj • sign ((— V/(x('^))) .) 



Let a :— ^ 



k+2 



Update := x^^^ + a{s — x^^^) 

end for 



Observe that in each iteration, this algorithm only introduces at most one new non-zero 
coordinate, so that the sparsity of x^'^^ is always upper bounded by the number of steps k. 
This means that the method is again of coordinate-descent-type, as in the simplex case of the 
previous Section 4.1. Its convergence analysis again directly follows from the general analysis 
from Section 3. 

Theorem 13 (Convergence of Sparse Greedy on the ^i-Ball). For each k > 1, the iterate x^^'> of 
Algorithm 4 satisfies 

where x* £ (}n is an optimal solution to problem (18). 

Furthermore, for any e > 0, after at most 2 + ~ ^ (e) ''^(^f^U steps, it has an iterate 

x^^^ of sparsity O (^), satisfying g{x^^'^) < e. 

Proof. This is a corollary of Theorem 3 and Theorem 5. □ 

The Duality Gap, and Duality of the Norms. We recall the definition of the duality gap (5) 
given by the linearization at any point x G On, see Section 2. Thanks to our Observation 2, the 
computation of the duality gap in the case of the £i-ball here becomes extremely simple, and is 
given by the norm that is dual to the £i-norm, namely the ^oo-norm of the used subgradient, 
i.e., 

g{x,dx) = \\dx\\oo +x^dx, and 
g(x) = ||V/(x)|L+x^V/(x) . 

Alternatively, the same expression can also be derived directly (without explicitly using duality 
of norms) by applying the Observation 12. 
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A Lower Bound on the Sparsity. The lower bound of (^) on the sparsity as proved in 
Section 4.2 for the simplex case in fact directly translates to the ^i-ball as well. Instead of 
choosing the objective function / as the distance to the origin (which is part of the £i-ball), 
we consider the optimization problem min /(x) := ||x — r\\2 with respect to the fixed point 

|jx||j<l 

r := (^,...,^) E M". This problem is of the form (18), and corresponds to optimizing the 
Euclidean distance to the point r given by mirroring the origin at the positive facet of the ii- 
ball. Here by the "positive facet" , we mean the hyperplane defined by the intersection of the 
boundary of the £i-ball with the positive orthant, which is exactly the unit simplex. Therefore, 
the proof for the simplex case from Section 4.2 holds analogously for our setting here. 

We have thus obtained that sparsity O (^) as obtained by the greedy Algorithm 4 is indeed 
best possible for £i-regularized optimization problems of the form (18). 



Using Barycentric Coordinates Instead. Clarkson [ClalO, Theorem 4.2] already observed that 
Algorithm 3 over the simplex A„ can be used to optimize a convex function f{y) over arbitrary 
convex hulls, by just using barycentric coordinates y = Ax, x E A„, for A E ]^"x»" being the 
matrix containing all m vertices of the convex domain as its columns. Here however we saw that 
for the ^i-ball, the steps of the algorithm are even slightly simpler, as well as that the duality 
gap can be computed instantly from the ^oo-norm of the gradient. 



Applications. Our Algorithm 4 applies to arbitrary convex vector optimization problems with 
an ||.||]^-norm regularization term, giving a guaranteed sparsity of O (^) for all these applications. 

A classical example for problems of this class is given by the important \\.\\-^^-regularized least 
squares regression approach, i.e. 

min \\Ax — b\\l + u \\x\\-, 

for a fixed matrix A E M"*^"", a vector b E M™ and a fixed regularization parameter fi > 
0. The same problem is also known as basis pursuit de-noising in the compressed sensing 
literature, which we will discuss more precisely in Section 5.1. The above formulation is in fact 
the Lagrangian formulation of the corresponding constrained problem for \\x\\i < t for some 
fixed parameter t corresponding to /i. This equivalent formulation is also known as the Lasso 
problem [Tib96] which is 

min \\Ax — b\\o 

S.t. < t . 

The above formulation is exactly a problem of our above form (18), namely 

min \\tAx — b\\2 , 

if we rescale the argument x =: tx so that \\x\\-^ < 1. 

Another important application for our result is logistic regression with ||.||]^-norm regulariza- 
tion, see e.g. [KKB07], which is also a convex optimization problem [Rcn05]. The reduction to 
an ^i-problem of our form (18) works exactly the same way as described here. 



Related Work. As we mentioned above, the optimization problem (18) — if / is the squared 
error of a linear function — is very well studied as the Lasso approach, see e.g. [Tib96] and the 
references therein. For general objective functions / of bounded curvature, the above interesting 
trade-off between sparsity and the approximation quality was already investigated by [SSSZIO], 
and also by our earlier paper [GJ09] for the analogous case of optimizing over the simplex. 
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[SSSZIO, Theorem 2.4] shows a sparse convergence analogous to our above Theorem 13, for 
the "forward greedy selection" algorithm on problem (18), but only for the case that / is 
differentiable. 



5.1 Relation to Matching Pursuit and Basis Pursuit in Compressed Sensing 

Both our sparse greedy Algorithm 3 for optimizing over the simplex and also Algorithm 4 for 
general £i-problems are very similar to the technique of matching pursuit, which is one of the 
most popular techniques in sparse recovery in the vector case [Tro04]. 

Suppose we want to recover a sparse signal vector x G from a noisy measurement vec- 
tor Ax = y £ M"^. For a given dictionary matrix A G M™-^", matching pursuit iteratively 
chooses the dictionary element Ai G M*" that has the highest inner product with the cur- 
rent residual, and therefore reduces the representation error /(x) = \\Ax — y\\2 by the largest 
amount. This choice of coordinate i = argmax^ Aj[Ax — y) exactly corresponds'^ to the choice 
of i := argmiUj (V/(x^^^)) in Algorithm 3. 

Another variant of matching pursuit, called orthogonal matching pursuit (OMP) [Tro04, 
TG07], includes an extra orthogonalization step, and is closer related to the coreset algorithms 
that optimize over the all existing set of non-zero coordinates before adding a new one, see 
e.g. [ClalO, Algorithm 8.2], or the analogous "fully corrective" variant of [SSSZIO]. If y = Ax, 
with X sparse and the columns of A sufficiently incoherent, then OMP recovers the sparsest 
representation for the given y [Tro04]. 

The paper [Zhall] recently proposed another algorithm that generalizes OMP, comes with a 
guarantee on correct sparse recovery, and also corresponds to "completely optimize within each 



(V/(xW)). 



as 



coreset". The method uses the same choice of the new coordinate i := argmax^ 
in our Algorithm 4. However the analysis of [Zhall] requires the not only bounded curvature as 
in our case, but also needs strong convexity of the objective function (which then also appears 
as a multiplicative factor in the number of iterations needed). Our Algorithm 4 as well as the 
earlier method by [Zha03] are simpler to implement, and have a lower complexity per iteration, 
as we do not need to optimize over several currently non-zero coordinates, but only change one 
coordinate by a fixed amount in each iteration. 

Our Algorithm 4 for general £i-regularized problems also applies to solving the so called basis 
pursuit problem [CDS98, FNW07] and [BV04, Section 6.5.4], which is mini,.gRn \\x\\^ s.t. Ax = y. 
Note that this is in fact just the constrained variant of the corresponding "robust" -regularized 
least squares regression problem 

which is the equivalent trade-off variant of our problem of the form (18). [FNW07] propose a 
traditional gradient descent technique for solving the above least squares problem, but do not 
give a convergence analysis. 

Solution path algorithms with approximation guarantees for related problems (obtaining solu- 
tions for all values of the tradeoff parameter /i) have been studied in [GJLIO, GJL12a, GJL12b], 
and the author's PhD thesis [Jagll], building on the same duality gap concept we introduced 
in Section 2. 

*The objective function f(x) := \\Ax — y Hj can be written as f{x) = {Ax—y)'^(Ax—y) = x^A^Ax—2y^Ax—y^y, 
so its gradient is \7f{x) = 2A^ Ax - 2A^y = 2A'^ [Ax - y) G R". 
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6 Optimization with Bounded ^oo-Morm 

Applying our above general optimization framework for the special case of the domain being the 
ll-llo^^-norm unit ball, we again obtain a very simple greedy algorithm. The running time will 
again correspond to O (^) many gradient evaluations. Formally, we consider problems of the 
form 

minimize f(x) 

xeR" (19) 

s.t. \\x\\^ < 1 . 

We denote the feasible set, i.e. the H-Hj^-norm unit ball, by □„ := {x G \ \\x\\^ < 1}. For 
this set, it will again be very simple to implement the internal primitive operation of optimizing a 
linear function over the same domain. The following crucial observation allows us to implement 
ExactLinear (c, Dn) in a very simple way. This can also alternatively be interpreted as the 
known fact that the dual-norm to the ^oo-norm is the £i-norm, which also explains why the greedy 
algorithm we will obtain here is very similar to the ^i-version from the previous Section 5. 

Observation 14. For any vector c G M", it holds that 

r T 

s G argmax y c 

where G is the sign-vector of c, defined by the sign of each individual coordinate, i.e. 
{s% = sign(Q) G {-1,1}. 

Using this observation for c = —dx in our general Algorithm 1, we directly obtain the following 
simple method for optimization over a box-domain as depicted in Algorithm 5. 

Algorithm 5 Sparse Greedy on the Cube 

Input: Convex function /, target accuracy e 
Output: e-approximate solution for problem (19) 
Set := 
for /c = . . . oo do 

Compute the sign- vector s of V/(a;(^)), such that 

Si = sign((-V/(j;W)).) , i = l..n 
Let a := ^2 

Update x^'^+i) := x^^) + a(s - xW) 
end for 

The convergence analysis again directly follows from the general analysis from Section 3. 
Theorem 15. For each k>l, the iterate x^^^ of Algorithm 5 satisfies 

f{x^^'^) - fix*) < . 



where x* G is an optimal solution to problem (19). 



Furthermore, for any e > 0, after at most 2 
x'^^'^ with g{x^^^) < e. 



£ 



+ 1 = 0(1) many steps, it has an iterate 



Proof. This is a corollary of Theorem 3 and Theorem 5. □ 



26 



The Duality Gap, and Duality of the Norms. We recall the definition of the duality gap (5) 
given by the linearization at any point x G Dn, see Section 2. Thanks to our Observation 2, the 
computation of the duality gap in the case of the ^oo-ball here becomes extremely simple, and 
is given by the norm that is dual to the ^oo-norm, namely the ^i-norm of the used subgradient, 
i.e., 

g{x,dx) = \\dx\\i+ x'^dx, and 
g(:E) = ||V/(x)||i + x^V/(x) . 

Alternatively, the same expression can also be derived directly (without explicitly using duality 
of norms) by applying the Observation 14. 

Sparsity and Compact Representations. The analogue of "sparsity" as in Sections 4 and 5 in 
the context of our Algorithm 5 means that we can describe the obtained approximate solution x 
as a convex combination of few (i.e. O(^) many) cube vertices. This does not imply that x has 
few non-zero coordinates, but that we have a compact representation given by only O(^) many 
binary n- vectors indicating the corresponding cube vertices, of which x is a convex combination. 

Applications. Any convex problem under coordinate-wise upper and lower constraints can be 
transformed to the form (19) by re-scaling the optimization argument. A specific interesting 
application was given by [MRU], who have demonstrated that integer linear programs can 
be relaxed to convex problems of the above form, such that the solutions coincide with high 
probability under some mild additional assumptions. 

Using Barycentric Coordinates Instead. Clarkson [ClalO, Theorem 4.2] already observed that 
Algorithm 3 over the simplex A„ can be used to optimize a convex function f{y) over arbitrary 
convex hulls, by just using barycentric coordinates y = Ax, x S A„, for A G ]R"x™- being the 
matrix containing all m vertices of the convex domain as its columns. Here however we saw 
that for the unit box, the steps of the algorithm are much simpler, as well as that the duality 
gap can be computed instantly, without having to explicitly deal with the exponentially many 
vertices (here m = 2") of the cube. 

7 Semidefinite Optimization with Bounded Trace 

We will now apply the greedy approach from the previous Section 3 to semidefinite optimiza- 
tion problems, for the case of bounded trace. The main paradigm in this section will be to 
understand the best achievable low-rank property of approximate solutions as a function of the 
approximation quality. 

In particular, we will show that our general Algorithm 1 and its analysis do lead to Kazan's 
method for convex semidefinite optimization with bounded trace, as given by [Haz08]. Kazan's 
algorithm can also be used as a simple solver for general SDPs. [Haz08] has already shown that 
guaranteed e-approximations of rank O (^) can always be found. Here we will also show that 
this is indeed optimal, by providing an asymptotically matching lower bound in Section 7.4. 
Furthermore, we fix some problems in the original analysis of [Haz08], and require only a weaker 
approximation quality for the internal linearized primitive problems. We also propose two 
improvement variants for the method in Section 7.3. 

Later in Section 11, we will discuss the application of these algorithms for nuclear norm and 
max-norm optimization problems, which have many important applications in practice, such as 
dimensionality reduction, low-rank recovery as well as matrix completion and factorizations. 
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We now consider convex optimization problems of the form (1) over the space X = g"^" of 
symmetric matrices, equipped with the standard Frobenius inner product {X, Y) = X . It is 

2 

left to the choice of the reader to identify the symmetric matrices either with M" and consider 
functions with f{X) = f{X'^), or only "using" the variables in the upper right (or lower left) 
triangle, corresponding to M"("+^)/2. In any case, the subgradients or gradients of our objective 
function / need to be available in the same representation (same choice of basis for the vector 
space X). 

Formally, we consider the following special case of the general optimization problems (1), i.e., 

minimize f(X) 

s.t. Tr(X) = 1 , (20) 
X 

We win write S := {X e S"""" | ^ ^ 0, Tr(X) = 1} for the feasible set, that is the PSD matrices 
of unit trace. The set S is sometimes called the Spectahedron, and can be seen as a natural 
generalization of the unit simplex to symmetric matrices. By the Cholesky factorization, it can 
be seen that the Spectahedron is the convex hull of all rank-1 matrices of unit trace (i.e. the 
matrices of the form vv'^ for a unit vector v E M", ||u||2 = 1). 

7.1 Low-Rank Semidefinite Optimization with Bounded Trace: The O(^) 
Algorithm by Hazan 

Applying our general greedy Algorithm 1 that we studied in Section 3 to the above semidef- 
inite optimization problem, we directly obtain the following Algorithm 6, which is Kazan's 
method [Haz08, GMll]. 

Note that this is now a first application of Algorithm 1 where the internal linearized problem 
ApproxLinear() is not trivial to solve, contrasting the applications for vector optimization 
problems we studied above. The algorithm here obtains low-rank solutions (sum of rank-1 
matrices) to any convex optimization problem of the form (20). More precisely, it guarantees e- 
small duality gap after at most O (^) iterations, where each iteration only involves the calculation 
of a single approximate eigenvector of a matrix M G S"^"-. We will see that in practice for 
example Lanczos' or the power method can be used as the internal optimizer ApproxLinear(). 

Algorithm 6 Kazan's Algorithm / Sparse Greedy for Bounded Trace 

Input: Convex function / with curvature constant Cf, target accuracy £ 
Output: e-approximate solution for problem (20) 
Set A:(o) := vv'^ for an arbitrary unit length vector v G M". 
for A: = . . . oo do 
Let a := ^ 

Compute V := v^^^ = ApproxEV {V f{X^^^),aCf) 
Update X^'^+i) := + a{vv^ - X^^)) 
end for 



Kere ApproxEV(^, e') is a subroutine that delivers an approximate smallest eigenvector (the 
eigenvector corresponding to the smallest eigenvalue) to a matrix A with the desired accuracy 
e' > 0. More precisely, it must return a unit length vector v such that Av < Amin(^) + e'- 
Note that as our convex function / takes a symmetric matrix X as an argument, its gradients 
V/(X) are given as symmetric matrices as well. 

If we want to understand this proposed Algorithm 6 as an instance of the general convex 
optimization Algorithm 1, we just need to explain why the largest eigenvector should indeed 
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be a solution to the internal linearized problem ApproxLinear(), as required in Algorithm 1. 
Formally, we have to show that v := ApproxEV(A, e') does approximate the linearized problem, 
that is 

vv'^ • A < min Y • A + e' 

Yes 

for the choice of v := ApproxEV(A, e'), and any matrix A G S"^". 

This fact is formalized in Lemma 16 below, and will be the crucial property enabling the fast 
implementation of Algorithm 6. 

Alternatively, if exact eigenvector computations are available, we can also implement the exact 
variant of Algorithm 1 using ExactLinear(), thereby halving the total number of iterations. 

Observe that an approximate eigenvector here is significantly easier to compute than a pro- 
jection onto the feasible set S. If we were to find the ||.||^^Q-closest PSD matrix to a given 
symmetric matrix A, we would have to compute a complete eigenvector decomposition of A, and 
only keeping those corresponding to positive eigenvalues, which is computationally expensive. 
By contrast, a single approximate smallest eigenvector computation as in APPR0XEV(A, e') 
can be done in near linear time in the number of non-zero entries of A. We will discuss the 
implementation of APPR0XEV(^, e') in more detail further below. 

Sparsity becomes Low Rank. As the rank-1 matrices are indeed the "vertices" of the domain S 
as shown in Lemma 16 below, our Algorithm 6 can be therefore seen as a matrix generalization 
of the sparse greedy approximation algorithm of [ClalO] for vectors in the unit simplex, see 
Section 4, which has seen many successful applications. Here sparsity just gets replaced by 
low rank. By the analysis of the general algorithm in Theorem 3, we already know that we 
obtain e-approximate solutions for any convex optimization problem (20) over the spectahedron 
S. Because each iterate X^^^ is represented as a sum (convex combination) of k many rank-1 
matrices vv'^ , it follows that X^'^^ is of rank at most k. Therefore, the resulting e- approximations 
are of low rank, i.e. rank 0(^)- 

For large-scale applications where ^ <^ n, the representation of X^*^^ as a sum of rank-1 
matrices is much more efficient than storing an entire matrix 

j^(fc) ^ gnxn_ Lg^^gj. Section 11.5 
(or see also [JSIO]) we will demonstrate that Algorithm 6 can readily be applied to practical 
problems for n > 10^ on an ordinary computer, well exceeding the possibilities of interior point 
methods. 

[Haz08] already observed that the same Algorithm 6 with a well-crafted function / can also be 
used to approximately solve arbitrary SDPs with bounded trace, which we will briefiy explain 
in Section 7.2. 

Linearization, the Duality Gap, and Duality of the Norms. Here we will prove that the general 
duality gap (5) can be calculated very efficiently for the domain being the spectahedron S. From 
the following Lemma 16, we obtain that 

g{X) =X . V/(X) + A^ax(-V/(X)) 

=X.V/(X)-A^in(V/(X)) . 

As predicted by our Observation 2 on formulating the duality gap, we have again obtained the 
dual norm to the norm that determines the domain D. It can be seen that over the space of 
symmetric matrices, the dual norm of the matrix trace-norm (also known as the nuclear norm) 
is given by the spectral norm, i.e. the largest eigenvalue. To see this, we refer the reader to the 
later Section 11.2 on the properties of the nuclear norm and its dual characterization. 

The following Lemma 16 shows that any linear function attains its minimum and maximum 
at a "vertex" of the Spectahedron S, as we have already proved for the case of general convex 
hulls in Lemma 8. 
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Lemma 16. The spectahedron is the convex hull of the rank-1 matrices, 

S = couv({vv'^ 1 G M", WvW^ = l}) • 

Furthermore, for any symmetric matrix A S S"^", it holds that 

maxA* X = Amax(^) • 
x&s 

Proof. Clearly, it holds that vv"^ G S for any unit length vector v £ M", as Tr:{vv'^) = \\v\\2. 
To prove the other inclusion, we consider an arbitrary matrix X £ S, and let X = U^U be its 
Cholesky factorization. We let ai be the squared norms of the rows of U, and let Ui be the row 
vectors of U, scaled to unit length. From the observation 1 = Tr(X) = Tr^U'^U) = Tr{UU'^) = 
- ai it follows that any X £ S can be written as a convex combination of at most n many 
rank-1 matrices X = X^ILi ^i'^i'^J with unit vectors Ui S M", proving the first part of the claim. 
Furthermore, this implies that we can write 

n n 

max A» X = max A • aiUiuf = max \^ ai{A • UiuJ), 

where the maximization max^.^^^ is taken over unit vectors Uj G M", ||n,|| = 1, for 1 < i < n, 
and real coefficients > 0, with Y^=\ Q^i = 1- Therefore 

max A* X = max a^iA • UjuT) 

i=l 

= max A • vv'^ 
veRi^ ,\\v\\=i 

= max Av 

t)eK",||i)||=i 

— -^max (^) ) 

where the last equality is the variational characterization of the largest eigenvalue. □ 



Curvature. We know that the constant in the actual running time for a given convex function 
/ : S'^^'* _). ]^ ig given by the curvature constant Cj as given in (7), which for the domain S 
becomes 

Cf := sup ^ (/(y) - fix) + iY-X). VfiX)) . (22) 
x,ye<s,ae[o,i], 

Y=X+a{V-X) 

Convergence. We can now see the convergence analysis for Algorithm 6 following directly as 
a corollary of our simple analysis of the general framework in Section 3. The following theorem 
proves that 0{l) many iterations are sufficient to obtain primal error < e. This result was 
already known in [Haz08, Theorem 1], or [GMll, Chapter 5] where some corrections to the 
original paper were made. 

Theorem 17. For each k >1, the iterate X^^^ of Algorithm 6 satisfies 

f{X^^))-f{X*)<j-^. 

where X* £ S is an optimal solution to problem (20). 

Furthermore, for any e > 0, after at most 2 + 1 = O (^) many steps, it has an iterate 

X^^"* of rank O [\) , satisfying g{X^^'^) < e. 

Proof. This is a corollary of Theorem 3 and Theorem 5. □ 
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Approximating the Largest Eigenvector. Approximating the smallest eigenvector of a sym- 
metric matrix V/(X) (which is the largest eigenvector of —Vf{X)) is a well-studied problem in 
the literature. We will see in the following that the internal procedure ApproxEV(M, e'), can be 
performed in near-linear time, when measured in the number of non-zero entries of the gradient 
matrix V/(X). This will follow from the analysis of [KW92] for the power method or Lanczos' 
algorithm, both with a random start vector. A similar statement has been used in [AHK05, 
Lemma 2]. 

Theorem 18. Let M £ S"X" &e a positive semidefinite matrix. Then with high probability, both 
i) O ^ '"^"^ ^ iterations of the power method, or 

j^j^j Q ^ ^ iterations of Lanczos' algorithm 

will produce a unit vector x such that ^^^^^ > 1 — 7. 

Proof. The statement for the power method follows from [KW92, Theorem 3.1(a)], and for 
Lanczos' algorithm by [KW92, Theorem 3.2(a)]. □ 



The only remaining obstacle to use this result for our internal procedure ApproxEV(-/Vf, e') is 
that our gradient matrix M = — V/(X) is usually not PSD. However, this can easily be fixed by 
adding a large enough constant t to the diagonal, i.e. M := M -\- tl, or in other words shifting 
the spectrum of M so that the eigenvalues satisfy Aj(M) = \i{M) + t > Q \/i. The choice of 
t = — \rain{M) \s good enough for this to hold. 

Now by setting 7 := < ^' ~ for some upper bound L > Amax(M) = Ainax(M) - Amin(^), 
this implies that our internal procedure ApproxEV(M, e') can be implemented by performing 
O many Lanczos steps (that is matrix- vector multiplications). Note that a simple 

choice for L is given by the spectral norm of M, since 2 ||M||gpg^ = 2maxj Aj(M) > Amax(-^) — 
Amin(-^)- We state the implication for our algorithm in the following corollary. 

Theorem 19. For M € S"^", and e' > 0, the procedure ApproxEV{M, e') requires a total of 
Q ^jy^ logWv^ ^ many arithmetic operations, with high probability, by using Lanczos' algorithm. 

Here Nf is the number of non-zero entries in M, which in the setting of Algorithm 6 is the 
gradient matrix — V/(Ar). We have also assumed that the spectral norm of M is bounded by L. 

Since we already know the number of necessary "outer" iterations of Algorithm 6, by Theo- 
rem 17, we conclude with the following analysis of the total running time. Here we again use 
that the required internal accuracy is given by e' = aCj < £Cf. 

Corollary 20. When using Lanczos ' algorithm for the approximate eigenvector procedure ApproxEV{., .), 
then Algorithm 6 provides an e-approximate solution in O (^) iterations, requiring a total of 

O ^^T3"^ arithmetic operations (with high probability). 

Here the notation 0{.) suppresses the logarithmic factor in n. This corollary improves the 
original analysis of [Haz08] by a factor of since [Haz08, Algorithm 1] as well as the proof 

of [Haz08, Theorem 1] used an internal accuracy bound of e' = O(j^) instead of the sufficient 
choice of e' = O(^) as in our general analysis here. 
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Representation of the Estimate X in the Algorithm. The above result on the total running 
time assumes the following: After having obtained an approximate eigenvector v, the rank-1 
update X^^^^'^ := (1 — a)X'^^'> + avv^ can be performed efficiently, or more precisely in time 
Nf. In the worst case, when a fully dense matrix X is needed, this update cost is Nf = n^. 
However, there are many interesting applications where the function / depends only on a small 
fraction of the entries of X, so that Nf <^ v?. Here, a prominent example is matrix completion 
for recommender systems. In this case, only those Nf many entries of X will be stored and 
affected by the rank-1 update, see also our Section 11.5. 

An alternative representation of X consists of the low-rank factorization, given by the v- 
vectors of each of the O(^) many update steps, using a smaller memory of size 0(7)- However, 
computing the gradient V/(X) from this representation of X might require more time then. 



7.2 Solving Arbitrary SDPs 

In [Haz08] it was established that Algorithm 6 can also be used to approximately solve arbitrary 
semidefinite programs (SDPs) in feasibility form, i.e., 

find X s.t. Ai • X < bi i = l..m 

XhO. ^ ' 

Also every classical SDP with a linear objective function 

maximize C • X 

X 

s.t. Ai»X<h i = l..m' (24) 
X^O . 

can be turned into a feasibility SDP (23) by "guessing" the optimal value C • X hy binary 
search [AHK05, Haz08]. 

Here we will therefore assume that we are given a feasibility SDP of the form (23) by its 
constraints Ai» X < bi, which we want to solve for X. We can represent the constraints of (23) 
in a smooth optimization objective instead, using the soft-max function 

/(X):=^log(^^e'^(^'-^-^')^ . (25) 

Suppose that the original SDP was feasible, then after O (^) many iterations of Algorithm 6, 
for a suitable choice of a, we have obtained X such that f{X) < e, which implies that all 
constraints are violated by at most e. This means that Ai • X < bi + e, or in other words we 
say that X is e-feasible [Haz08, GMll]. It turns out the best choice for the parameter a is 
'"^"^ , and the curvature constant Cf{a) for this function is bounded by a • maxj Amax(^i)^- The 
total number of necessary approximate eigenvector computations is therefore in O ^ ^"^^"^ ^ . In 
fact. Algorithm 6 when applied to the function (25) is very similar to the multiplicative weights 
method [AHK05]. Note that the soft-max function (25) is convex in X, see also [Ren05]. For a 
slightly more detailed exhibition of this approach of using Algorithm 6 to approximately solving 
SDPs, we refer the reader to the book of [GMll]. 

Note that this technique of introducing the soft-max function is closely related to smoothing 
techniques in the optimization literature [Nem04, BB09], where the soft-max function is intro- 
duced to get a smooth approximation to the largest eigenvalue function. The transformation to 
a smooth saddle-point problem suggested by [BB09] is more complicated than the simple notion 
of e-feasibility suggested here, and will lead to a comparable computational complexity in total. 
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7.3 Two Improved Variants of Algorithm 6 



Choosing the Optimal a by Line-Search. As we mentioned already for the general algorithm 
for convex optimization in Section 3, the optimal a in Algorithm 6, i.e. the a € [0, 1] of best 
improvement in the objective function / can be found by line-search. 

In particular for matrix completion problems, which we will discuss in more details in Sec- 
tion 11.5, the widely used squared error is easy to analyze in this respect: If the optimization 
function is given by /(X) = \ YlijeP^-^''i-i ~ ^u)^' where P is the set of observed positions of the 
matrix y, then the optimality condition (10) from Section 3.3 is equivalent to 

a = — 72 . 26 

Here X = X^^\ and v is the approximate eigenvector v^'^^ used in step k of Algorithm 6. The 
above expression is computable very efficiently compared to the eigenvector approximation task. 



Immediate Feedback in the Power Method. As a second improvement, we propose a heuristic 
to speed up the eigenvector computation, i.e. the internal procedure ApproxEV (V/(X), e'). 
Instead of multiplying the current candidate vector with the gradient matrix V/(X) in each 
power iteration, we multiply with ^ -|- V/(X)), or in other words the average between 

the current gradient and the gradient at the new candidate location X = {l-l)X^''^ + lv^''^v'^''^ . 
Therefore, we immediately take into account the effect of the new feature vector v^'^^ This 
heuristic (which unfortunately does not fall into our current theoretical guarantee) is inspired by 
stochastic gradient descent as in Simon Funk's method, which we will describe in Section 11.5.4. 
In practical experiments, this proposed slight modification will result in a significant speed-up 
of Algorithm 6, as we will observe e.g. for matrix completion problems in Section 11.5. 



7.4 r2(^) Lower Bound on the Rank 

Analogous to the vector case discussed in Section 4.2, we can also show that the rank of O(^), 
as obtained by the greedy Algorithm 6 is indeed optimal, by providing a lower bound of f^(^). 
In other words we can now exactly characterize the trade-off between rank and approximation 
quality, for convex optimization over the spectahedron. 

For the lower bound construction, we consider the convex function f{X) := ||X||p'ro = X • X 
over the symmetric matrices S"^". This function has gradient V/(X) = 2X. We will later see 
that its curvature constant is Cj = 2. 

The following lemmata show that the above sparse SDP Algorithm 6 is optimal for the ap- 
proximation quality (primal as well as dual error respectively), giving lowest possible rank, up 
to a small multiplicative constant. 

Lemma 21. For f{X) := \\X\\'^p^^, and 1 < k < n, it holds that 

min f(X) = —. 
Rk{X)<k 

We will see that this claim can be reduced to the analogous Lemma 10 for the vector case, 
by the standard technique of diagonalizing a symmetric matrix. (This idea was suggested by 
Elad Hazan) . Alternatively, an explicit (but slightly longer) proof without requiring the spectral 
theorem can be obtained by using the Cholesky-decomposition together with induction on k. 
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1 1 2 

Proof. We observe that the objective function H-H^^o, the trace, as well as the property of being 
positive semidefinite, are all invariant under orthogonal transformations (or in other words under 
the choice of basis). 

By the standard spectral theorem, for any symmetric matrix X of Rk(X) < k, there exists an 
orthogonal transformation mapping X to a diagonal matrix X' with at most k non-zero entries 
on the diagonal (being eigenvalues of X by the way). For diagonal matrices, the IMI^j.^ matrix 
norm coincides with the ||.||2 vector norm of the diagonal of the matrix. Finally by applying the 
vector case Lemma 10 for the diagonal of X', we obtain that f{X) = f{X') > p 

To see that the minimum can indeed be attained, one again chooses the "uniform" example 
X := ^Ifc G S, being the matrix consisting of k non-zero entries (of ^ each) on the diagonal. 
This gives f{X) = l. □ 

Recall from Section 7.1 that for convex problems of the form (20) over the Spectahedron, the 
duality gap is the non-negative value g{X) := f{X) — uj{X) = X»V/(X) — Amin(V/(X)). Also, 
by weak duality as given in Lemma 1, this value is always an upper bound for the primal error, 
that is f{X) - f{X*) < g{X) VX. 

Lemma 22. For f{X) := \\X\\^p^^, and any k £ N, k < n, it holds that 

g{X) >\ \fX eS s.t. Rk{X) < k. 

Proof. g{X) = Amax(-V/(X))+X.V/(X) = -Amin(X)+X.2X. Wenowusethat Amm(X) = 
for all symmetric PSD matrices X that are not of full rank n, and that by Lemma 21, we have 
XmX = Ti{X'^X)= f[X)>\. □ 

The Curvature. We will compute the curvature Cf of our function f{X) := X • X, showing 
that Cj = 2 in this case. Using the definition (7), and the fact that here 

f{Y)-f{X)-{Y-X).Vf{X) 
= Y •¥ - X» X - {Y - X)»2X 

= \\X — Y\\p^^ , 

one obtains that Cf = sup^ ||X — = (i\a,m.pro{SY = 2. Finally the following 

Lemma 23 shows that the diameter is indeed 2. 

Lemma 23 (Diameter of the Spectahedron). 

dia.m.FroiS)'^ = 2 . 

Proof. Using the fact that the spectahedron S is the convex hull of the rank-1 matrices of unit 

trace, see Lemma 16, we know that the diameter must be attained at two vertices of S, i.e. 

u,v £ W" with ||n||2 = ||u||2 = Ij and 

II T Tl|2 
\\vv — uu „ 

II II Fro 

= ||„||4+ ||^||4_2(^T^)2 , 

Clearly, this quantity is maximized if u and v are orthogonal. □ 

Note: We could also study f{X) := 7||X||^^^ instead, for some 7 > 0. This function has 
curvature constant Cf = 27, and for this scaling our above lower bounds will also just scale 
linearly, giving instead of ^. 
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8 Semidefinite Optimization with ^oo-Bounded Diagonal 



Here we specialize our general Algorithm 1 to semidefinite optimization problems where all 
diagonal entries are individually constrained. This will result in a new optimization method 
that can also be applied to max- norm optimization problems, which we will discuss in more 
detail in Section 11. As in the previous Section 7, here we also consider matrix optimization 
problems over the space X = g*^^" of symmetric matrices, equipped with the standard Frobenius 
inner product {X, Y) = X •¥ . 

Formally, we consider the following special case of the general optimization problems (1), i.e. 

minimize f(X) 

s.t. Xii < 1 Vi, (27) 

X to . 

We will write S := {X £ S"^" \ X tO, Xu < 1 Vi} for the feasible set in this case, that is the 
PSD matrices whose diagonal entries are all upper bounded by one. This class of optimization 
problems has become widely known for the linear objective case when f{X) = A»X, if A being 
the Laplacian matrix of a graph. In this case, one obtains the standard SDP relaxation of the 
Max-Cut problem [GW95], which we will briefly discuss below. Also, this optimization domain 
is strongly related to the matrix max-norm, which we study in more detail in Section 11.3. 

Our general optimization Algorithm 1 directly applies to this specialized class of optimization 
problems as well, in which case it becomes the method depicted in the following Algorithm 7. 

Algorithm 7 Sparse Greedy for Max-Norm Bounded Semidefinite Optimization 
Input: Convex function / with curvature Cf, target accuracy e 
Output: e-approximate solution for problem (27) 
Set := vv'^ for an arbitrary unit length vector v G M". 
for A; = . . . oo do 
Let a := ^2 

Compute S := ApproxLinear (V/(A:W), ffl, aC7j) 
Update X^'^+i) := AT^ + a{S - X^^^) 
end for 



The Linearized Problem. Here, the internal subproblem ApproxLinear() of approximately 
minimizing a linear function over the domain ffl of PSD matrices is a non-trivial task. Every 
call of ApproxLinear(A, ffl, e') in fact means that we have to solve a semidefinite program 
minyga Y • A for a given matrix yl, or in other words 

minimize Y • A 

Y 

S.t. Yii < 1 Vi, (28) 
Y to 

up to an additive approximation error of e' = aC j . 

Relation to Max-Cut. In [AHK05, Kal07], the same linear problem is denoted by (MaxQP). 
In the special case that A is chosen as the Laplacian matrix of a graph, then the above SDP 
is widely known as the standard SDP relaxation of the Max-Cut problem [GW95] (not to be 
confused with the combinatorial Max-Cut problem itself, which is known to be NP-hard). In 
fact the original relaxation uses equality constraints Ya = 1 on the diagonal instead, but for any 
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matrix A of positive diagonal entries (such as e.g. a graph Laplacian), this condition follows 
automatically in the maximization variant of (28), see [KL96], or also [GMll, Kal07] for more 
background. 

Duality and Duality of Norms. In Section 11.3 we will see that the above quantity (28) that 
determines both the step in our greedy Algorithm 7, but also the duality gap, is in fact the norm 
of A that is dual to the matrix max-norm. 

For optimization problems of the form (27), it can again be shown that the poor-man's duality 
given by the linearization (see also Section 2) indeed coincides with classical Wolfe- duality from 
the optimization literature. 

Fortunately, it was shown by [AHK05] that also this linearized convex optimization prob- 
lem (28) — and therefore also our internal procedure ApproxLinear(.) — can be solved rela- 
tively efficiently, if the matrix A (i.e. V f{X) in our case) is sparse.^ 

Theorem 24. The algorithm of [AHK05] delivers an additive e' -approximation to the linearized 
problem (28) in time 

where the constant L > is an upper bound on the maximum value ofY»A over y € ffl, and 
Na is the number of non-zeros in A. 

Proof. The results of [AHK05, Theorem 3] and [Kal07, Theorem 33] give a running time of order 
O • min |a^, ^° obtain a multiplicative (1 — e)-approximation, where a* is the value 

of an optimal solution. Formally we obtain S" G ffl with S • A > (1 — e)a* . In other words by 
using an accuracy of e := we obtain an additive e'-approximation to (28). □ 

Here the notation 0(.) again suppresses poly-logarithmic factors in n, and N is the number 
of non-zero entries of the matrix A. Note that analogous to the approximate eigenvector com- 
putation for Hazan's Algorithm 6, we need the assumption that the linear function given by 
y • V/(X) is bounded over the domain y G ffl. However this is a reasonable assumption, as our 
function has bounded curvature Cf (corresponding to V/(X) being Lipschitz-continuous over 
the domain ffl), and the diameter of ffl is bounded. 

The reason we need an absolute approximation quality lies in the analysis of Algorithm 1, 
even if it would feel much more natural to work with relative approximation quality in many 
cases. 

Convergence. The convergence result for the general Algorithm 1 directly gives us the analysis 
for the specialized algorithm here. Note that the curvature over the domain ffl here is given by 

Cf := sup ^ (/(y) - f{X) + {Y-X). Vf{X)) . (29) 

x,yeffl,ae[o,i], 

Y=X+aiV-X) 

Theorem 25. For each k >1, the iterate X^'^^ of Algorithm 7 satisfies 

f{X('))-f{X*)<^. 



^Also, Kale in [Kal07, Theorem 14] has shown that this problem can be solved very efficiently if the matrix 
A — — V/(X) is sparse. Namely if A is the Laplacian matrix of a weighted graph, then a multiplicative 
£-approximation to (28) can be computed in time 0{^Na) time, where Na is the number of non-zero entries 
of the matrix A. Here A is the maximum entry on the diagonal of A, and d is the average value on the 
diagonal. 
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where X* £ S is an optimal solution to problem (27). 



Furthermore, after at most 2 
< e. 



8Cf 



+ 1 



O(-) many steps, it has an iterate X^'^) with 



Proof. This is a corohary of Theorem 3 and Theorem 5. 



□ 



Applications. The new algorithm can be used to solve arbitrary max-norm constrained convex 
optimization problems, such as max-norm regularized matrix completion problems, which we 
will study in Section 11. 



9 Sparse Semidefinite Optimization 

Another interesting optimization domain among the semidefinite matrices is given by the ma- 
trices with only one non-zero off-diagonal entry. Here we specialize our general Algorithm 1 to 
convex optimization over the convex hull given by such matrices. Our algorithm will therefore 
obtain e-approximate solutions given by only O (^) such sparse matrices, or in other words 
solutions of sparsity O (^)- 

Why bother? The same sparse matrices are also used in the graph sparsification approach 
by [BSS09]6. Furthermore, sparse solutions to convex matrix optimization problems have gained 
interest in dimensionality reduction, as in sparse PCA, see [ZdGlO] for an overview. 



Setup. Formally, here we again use the standard Frobenius inner product {X, Y) = X • Y 
over the symmetric matrices S"^", and consider the sparse PSD matrices given by Pfe) := 

for any fixed pair of indices i,j £ [n], i ^ j- In other words 



(ei + ej){ei + e^)^ = 

Puv^ = 1 for u G {i,j},v G and zero everywhere else. We also consider the analogous 



"negative" counterparts of such matrices, namely N^'''"'' := (e^ — ej)(ej — 



1 

-1 



-1 
1 



i.e. NuJ^ = —1 for the two off-diagonal entries {u,v) £ {{i,j), ij,i)}, and NuJ^ = 1 for the two 
diagonal entries (u, f) G {{i,i), and zero everywhere else. 

Analogously to the two previous applications of our method to semidefinite optimization, we 
now optimize a convex function, i.e. 



minimize f{X) (30) 

over the domain 



D = 5,Vse := conv [J U [J N^'^^ 



Optimizing over Sparse Matrices, and Solving the Linearization. Applying our general Al- 
gorithm 1 to this class of problems (30) becomes very simple, as the linear primitive problem 
ExactLinear (-Dx ) ^spjjj.gg) for any fixed matrix Dx G g"^*^ is easily solved over 3^^^^^^^,. From 
our simple Lemma 8 on linear functions over convex hulls, we know that this linear minimum 
is attained by the single sparse matrix 

p{ij) or iV(^^) that maximizes the inner product with 



^The theoretical result of [BSS09] guarantees that all eigenvalues of the resulting sparse matrix (corresponding 
to the Laplacian of a sparse graph) do not differ too much from their counterparts in the original graph. 
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—Dx- The optimal pair of indices {k,l) can be found by a linear pass through the gradient 
Dx = V/(X). This means that the linearized problem is much easier to solve than in the above 
two Sections 7 and 8. Altogether, Algorithm 1 will build approximate solutions X'^^^ each of 
which is a convex combination of k of the atomic matrices P^*-') or N^'^^\ as formalized in the 
following theorem: 

Theorem 26. Let n > 2 and let X^^^ := P^^^) be the starting point. Then for each k > I, the 
iterate X^'^^ of Algorithm 1 has at most 4(fc + 1) non-zero entries, and satisfies 

f{X^^))-f{X*)<j-^. 

where X* G 5^parse ^-^ optimal solution to problem (30). 

Furthermore, for any e > Q, after at most 2 + 1 = O (^) many steps, it has an iterate 

X(^) of only O i^) many non-zero entries, satisfying g{X^^^) < e. 

Proof. This is a corollary of Theorem 3 and Theorem 5. The sparsity claim follows from our 
observation that the step directions given by ExactLinear (V/(Ar), 5gp3^j.gg) are always given 
by one of the sparse matrices P^*-') or N^^^\ □ 



Optimizing over Non-Negative Matrix Factorizations. We also consider the slight variant 
of (30), namely optimizing only over one of the two types of matrices as the domain D, i.e. only 
combinations of positive p(*-j) or only of negative A^fe'). This means that the domain is given 
hy D = 5^ptij,gg := conv (^Ujj -f''*"'^) or D = 5^p^j,gg := conv (^Ujj • The above analysis for 

Algorithm 1 holds in exactly the same way. Now for S^^^^.^^, each step direction s = s^*^) used 
by the algorithm is given by s = P(^J) = (Bi + ej){ei + ejf for some i, 7, and so we have that 
each of the approximations X^*^) is a sum of k many positive rank-1 factors of this form. In 
other words in each step k, X^^') = LR^ is a product of two (entry-wise) non-negative matrices 
of at most k columns each, i.e. L G R"^'^ and P"^'^. Consequently, our algorithm provides 
solutions that are non-negative matrix factorizations, which is a successful technique in matrix 
completion problems from recommender systems, see e.g. [Wu07]. 



Relation to Bounded Trace and Diagonally Dominant Matrices. Observe the matrices in 
'^^'arse ^ovui a subset of the bounded trace PSD matrices S that we studied in the previous 
Section 7, since every matrix p(*-j) or A^^*-?) is PSD and has trace equal two. Furthermore, we 
observe that all matrices X G ^^parse ^-re diagonally dominant, meaning that 

\Xii\>Y,\Xij\ yie[n] 

In the case that we restrict to using only one of the two types of matrices S^p^^^.^^ or S^~^j.g^ as 
the domain, then we have that equality |Ajj| = X^^yj always holds, since this equality is 
preserved under taking convex combinations, and holds for the atomic matrices P^*-') and N^'^^\ 



The Curvature. The above reasoning also implies that the curvature Cj for problems of the 
form (30) is upper bounded by the curvature in the spectahedron-case as given in (22), since 

c4+ c 9^ C 9 • 9 
'-^sparse — '-^sparse — ^ • 
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Applications and Future Research. Computationally, the approach here looks very attractive, 
as the cost of a "sparse" step here is much cheaper than an approximate eigenvector computation 
which is needed in the bounded trace explained in Section 7. 

Also, it will be interesting to see how a regularization by constraining to a scaled domain 
'^sparse '^^parse ^^^^ perform in practical machine learning applications as for example dimen- 
sionality reduction, compared to nuclear norm regularization that we will discuss in the following 
Chapter 11. 

It also remains to investigate further on whether we can approximate general bounded trace 
semidefinite problems of the form (20) by using only sparse matrices. 



10 Submodular Optimization 

For a finite ground set S, a real valued function defined on all subsets of S, is called submodular, 
if 

g{XnY)+g{XUY)<g{X)+g{Y) VX,yc5 

For any given submodular function g with g(f/>) = 0, the paper [Lov83, Section 3] introduces a 
corresponding convex set in M}^^, called the submodular polyhedron (or also Lovasz polyhedron). 



We would now like to study convex optimization problems over such domains, which become 
compact convex sets if we intersect with the positive orthant, i.e. D := VgD M>o. 

Nicely for our optimization framework, [Lov83, Section 3] already showed that there is a 

simple greedy algorithm which optimizes any linear function over the domain Vg, i.e. it solves 

maxc^x, or in other words it exactly solves our internal problem ExactLinear (c, Pg). 
xeVg 

Lovasz [Lov83] already demonstrated how to use this kind of linear optimization over Vg to 
solve submodular minimization problems. It remains to investigate if there are interesting ap- 
plications for the wider class of more general convex (non-linear) functions / over such domains, 
as addressed by our Algorithm 1. 



11 Optimization with the Nuclear and Max-Norm 

Matrix optimization problems with a nuclear norm or max-norm regularization, such as e.g. 
low norm matrix factorizations, have seen many applications recently, ranging from low-rank 
recovery, dimensionality reduction, to recommender systems. We propose two new first-order 
approximation methods building upon two of the simple semidefinite optimizers we have studied 
above, that is the approximate SDP solver of [HazOS] from Section 7 on one hand, and our 
bounded diagonal optimizer from Section 8 on the other hand. The algorithms come with 
strong convergence guarantees. 

In contrast to existing methods, our nuclear norm optimizer does not need any Cholesky 
or singular value decompositions internally, and provides guaranteed approximations that are 
simultaneously of low rank. The method is free of tuning parameters, and easy to parallelize. 

11.1 Introduction 

Here we consider convex optimization problems over matrices, which come with a regularization 
on either the nuclear norm or the max-norm of the optimization variable. 
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Convex optimization with the nuclear norm has become a very successful technique in vari- 
ous machine learning, computer vision and compressed sensing areas such as low-rank recovery 
[FHBOl, CR09, CTIO], dimensionality reduction (such as robust principal component analy- 
sis [CLMWll]), and also recommender systems and matrix completion. Here matrix factor- 
izations [SRJ04, KBV09] — regularized by the nuclear or max-norm — have gained a lot of 
attention with the recently ended Netflix Prize competition. Many more applications of simi- 
lar optimization problems can be found among dimensionality reduction, matrix classification, 
multi-task learning, spectral clustering and others. The success of these methods is fueled by 
the property of the nuclear norm being a natural convex relaxation of the rank, allowing the use 
of scalable convex optimization techniques. 

Based on the semidefinite optimization methods that we have presented in the above Sections 7 
and 8, we propose two new, yet simple, first-order algorithms for nuclear norm as well as max- 
norm regularized convex optimization. 

For the nuclear norm case, our proposed method builds upon the first-order scheme for semidef- 
inite optimization by [Haz08], which we have investigated in Section 7.1. This approach allows us 
to significantly reduce the computational complexity per iteration, and therefore scale to much 
larger datasets: While existing methods need an entire and exact singular value decomposition 
(SVD) in each step, our method only uses a single approximate eigenvector computation per 
iteration, which can be done by e.g. the power method. A conference version of our work for 
nuclear norm regularized problems has appeared in [JSIO]. 

In the same spirit, we will also give a new algorithm with a convergence guarantee for optimiz- 
ing with a max-norm regularization. For matrix completion problems, experiments show that 
the max-norm can result in an improved generalization performance compared to the nuclear 
norm in some cases [SRJ04, LRS^IO]. 

Nuclear Norm Regularized Convex Optimization. We consider the following convex optimiza- 
tion problems over matrices: 

min f{Z)+fi\\Z\l (31) 



and the corresponding constrained variant 



min ^ fiZ) (32) 



mxrt ll^ll < 



where f{Z) is any diff'erentiable convex function (usually called the loss function), ||.||^ is the 
nuclear norm of a matrix, also known as the trace norm (sum of the singular values, or ^i-norm 
of the spectrum). Here fi > and t > respectively are given parameters, usually called the 
regularization parameter. 

The nuclear norm is know as the natural generalization of the (sparsity inducing) £i-norm 
for vectors, to the case of semidefinite matrices. When choosing f{X) := \\A{X) — b\\2 for some 
linear map A : M"-^™ — ). Rf, the above formulation (31) is the matrix generalization of the 
problem min^^gRn \\Ax — b\\2 + fJ-WxW-^, for a fixed matrix A, which is the important ii- regularized 
least squares problem, also known as the basis pursuit de-noising problem in the compressed 
sensing literature, see also Section 5.1. The analoguous vector variant of (32) is the Lasso 
problem [Tib96] which is min^igM" — b\\2 \\x\\i < t|. 

Max-Norm Regularized Convex Optimization. Intuitively, one can think of the matrix max- 
norm as the generalization of the vector ioQ-noim to PSD matrices. Here we consider optimiza- 
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tion problems with a max-norm regularization, which are given by 




max 



(33) 



and the corresponding constrained variant being 



mm 



fiz) . 



(34) 



max — 



Our Contribution. Using the optimization methods from the previous Sections 7 and 27, we 
present a much simpler algorithm to solve problems of the form (32), which does not need 
any internal SVD computations. The same approach will also solve the max-norm regularized 
problems (34). We achieve this by transforming the problems to the convex optimization setting 
over positive semidefinite matrices which we have studied in the above Sections 7.1 and 8. 

Our new approach has several advantages for nuclear norm optimization when compared 
to the existing algorithms such as "proximal gradient" methods (APG) and "singular value 
thresholding" (SVT), see e.g. [GLW+09, CCSIO, TYIO, JY09], and also in comparison to the 
alternating-gradient-descent-type methods (as e.g. [RS05, Lin07]). 

i) By employing the approximate SDP solver by [Haz08], see Algorithm 6, we obtain a 
guaranteed e-approximate solution Z after O(^) iterations. Crucially, the resulting solu- 
tion Z is simultaneously of low rank, namely rank OQ) . Also the algorithm maintains a 
compact representation of Z in terms of a low-rank matrix factorization Z = LR^ (with 
the desired bounded nuclear norm), and can therefore even be applied if the full matrix 
Z would be far too large to even be stored. 

ii) Compared to the alternating-gradient-descent-type methods from machine learning, we 
overcome the problem of working with non-convex formulations of the form /(LR^), 
which is NP-hard, and instead solve the original convex problem in f{Z). 

iii) The total running time of our algorithm for nuclear norm problems grows linear in the 
problem size, allows to take full advantage of sparse problems such as e.g. for matrix 
completion. More precisely, the algorithm runs in time O (^^^t^ j where Nf is the number 
of matrix entries on which the objective function / depends. Per iteration, our method 
consists of only a single approximate (largest) eigenvector computation, allowing it to 
scale to any problem size where the power method (or Lanczos' algorithm) can still be 
applied. This also makes the method easy to implement and to parallelize. Existing 
APG/SVT methods by contrast need an entire SVD in each step, which is significantly 
more expensive. 

iv) On the theory side, our simple convergence guarantee of O(^) steps holds even if the 
used eigenvectors are only approximate. In comparison, those existing methods that 
come with a convergence guarantee do require an exact SVD in each iteration, which 
might not always be a realistic assumption in practice. 

We demonstrate that our new algorithm on standard datasets improves over the state of the 
art methods, and scales to large problems such as matrix factorizations on the Netflix dataset. 

Kazan's Algorithm 6 can be interpreted as the generalization of the coreset approach to 
problems on symmetric matrices, which we have explained in the previous Section 7.1. Compared 
to the O^l/y/e) convergence methods in the spirit of [Ngs83, Nes07a], our number of steps is 
larger, which is however more than compensated by the improved step complexity, being lower 
by a factor of roughly (n -\- m) . 

Our new method for the nuclear norm case can also be interpreted as a modified, theoret- 
ically justified variant of Simon Funk's popular SVD heuristic [Web06] for regularized matrix 
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factorization. To our knowledge this is the first guaranteed convergence result for this class of 
alternating-gradient-descent-type algorithms. 

Related Work. For nuclear norm optimization, there are two lines of existing methods. On 
the one hand, in the optimization community, [TYIO, LST09], [GLW^09] and [JY09] indepen- 
dently proposed algorithms that obtain an e-accurate solution to (31) in 0{l/y/£) steps, by 
improving the algorithm of [CCSIO]. These methods are known under the names "accelerated 
proximal gradient" (APG) and "singular value thresholding" (SVT). More recently also [MHTIO] 
and [MGC09] proposed algorithms along the same idea. Each step of all those algorithms re- 
quires the computation of the singular value decomposition (SVD) of a matrix of the same size 
as the solution matrix, which is expensive even with the currently available fast methods such 
as PROPACK. [TYIO] and [JY09] and also [GLW+09] show that the primal error of their algo- 
rithm is smaller than e after 0{l/^/e) steps, using an analysis inspired by [Nes83] and [BT09]. 
For an overview of related algorithms, we also refer the reader to [CLMWll]. As mentioned 
above, the method presented here has a significantly lower computational cost per iteration (one 
approximate eigenvector compared to a full exact SVD), and is also faster in practice on large 
matrix completion problems. 

On the other hand, in the machine learning community, research originated from matrix 
completion and factorization [SRJ04] , later motivated by the Netflix prize challenge, getting sig- 
nificant momentum from the famous blog post by [Wcb06] . Only very recently an understanding 
has formed that many of these methods can indeed by seen as optimizing with regularization 
term closely related to the nuclear norm, see Section 11.5.4 and [JSIO, SSIO]. The majority of 
the currently existing machine learning methods such as for example [RS05, Lin07] and later 
also [Pat07, ZWSP08, KBV09, TPNT09, IRIO, GNHSll] are of the type of "alternating" gra- 
dient descent applied to f{LR^), where at each step one of the factors L and R is kept fixed, 
and the other factor is updated by a gradient or stochastic gradient step. Therefore, despite 
working well in many practical applications, all these mentioned methods can get stuck in local 
minima — and so are theoretically not well justified, see also the discussion in [DeC06] and our 
Section 11.4. 

The same issue also comes up for max-norm optimization, where for example [LRS^IO] op- 
timize over the non-convex factorization (38) for bounded max-norm. To our knowledge, no 
algorithm with a convergence guarantee was known so far. 

Furthermore, optimizing with a rank constraint was recently shown to be NP-hard [GGIO]. 
In practical applications, nearly all approaches for large scale problems are working over a 
factorization Z = LB?- of bounded rank, therefore ruling out their ability to obtain a solution 
in polynomial time in the worst-case, unless P = NP. 

Our new method for both nuclear and max-norm avoids all the above described problems by 
solving an equivalent convex optimization problem, and provably runs in near linear time in the 
nuclear norm case. 

11.2 The Nuclear Norm for Matrices 

The nuclear norm \\Z\\^ of a rectangular matrix Z G M'"^", also known as the trace norm or Ky 
Fan norm, is given by the sum of the singular values of Z, which is equal to the £i-norm of the 
singular values of Z (because singular values are always non- negative) . Therefore, the nuclear 
norm is often called the Schatten £i-norm. In this sense, it is a natural generalization of the 
£i-norm for vectors which we have studied earlier. 

The nuclear norm has a nice equivalent characterization in terms of matrix factorizations of Z, 
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I.e. 

||Z||, :=^min^-(||L||^^„+ , (35) 

where the number of columns of the factors L G ]^"»xfc and R^^^ is not constrained [FHBOl, 
SRJ04]. In other words, the nuclear norm constrains the average Euclidean row or column norms 
of any factorization of the original matrix Z. 

Furthermore, the nuclear norm is dual to the standard spectral matrix norm (i.e. the matrix 
operator norm), meaning that 

IIZIL = max B • Z , 

B,\\B\\s,ec<^ 

see also [RFPIO]. Recall that ll-Bll^pe^ is defined as the first singular value (7i{B) of the matrix B. 

Similarly to the property of the vector ||.||^-norm being the best convex approximation to the 
sparsity of a vector, as we discussed in Section 5 the nuclear norm is the best convex approx- 
imation of the matrix rank. More precisely, ||.||^ is the convex envelope of the rank [FHBOl], 
meaning that it is the largest convex function that is upper bounded by the rank on the convex 
domain of matrices |z 1 1-^1 1 spec — '^^i^ motivates why the nuclear norm is widely used 
as a proxy function (or convex relaxation) for rank minimization, which otherwise is a hard 
combinatorial problem. 

Its relation to semidefinite optimization — which explains why the nuclear norm is often called 
the trace norm — is that 

IIZIL = minimize t 

II* 

s,. (J, ^) to and (36) 
Tv{V) + Tr(iy) < 2t . 

Here the two optimization variables range over the symmetric matrices V S S'"^'" and W G 
gnxn^ This semidefinite characterization will in fact be the central tool for our algorithmic 
approach for nuclear norm regularized problems in the following. The equivalence of the above 
characterization to the earlier "factorization" formulation (35) is a consequence of the following 
simple Lemma 27. The Lemma gives a correspondence between the (rectangular) matrices 
Z G j^mxn q£ bounded nuclear norm on one hand, and the (symmetric) PSD matrices X G 
g(m+n)x(m+n) q£ bounded tracc on the other hand. 

Lemma 27 ([FHBOl, Lemma 1]). For any non-zero matrix Z G ]^rnxn t £M., it holds that 

\\z\l<- 

II II* - 2 



if and only if 



3 symmetric matrices V G S"^^"^, W G S'^^" 
s.t. w)-^ "^'^ '^^^^^ + Ti{W) < t . 



Proof. =^ Using the characterization (35) of the nuclear norm = min^^^r^^ ^(ll-^llFro + 

ll^llFro) we get that 3 L,R, LR^ = Z s.t. + \\Rfp^^ = TY(LL^) + Ti{RR^) < t, or in 

other words we have found a matrix ( ^ ^^j, ) = (^)(^)"^ ^ of trace < t. 
I <= I As the matrix ( ^ ) is symmetric and PSD, it can be (Cholesky) factorized to (L; R){L; R)^ 
s.t. LR^ = Z and t > Tr(LL^) + Tt{RR^) = ||L||^^„ + therefore \\Z\\^ < |. □ 
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Interestingly, for characterizing bounded nuclear norm matrices, it does not make any differ- 
ence whether we enforce an equality or inequality constraint on the trace. This fact will turn 
out to be useful in order to apply our Algorithm 6 later on. 

Corollary 28. For any non-zero matrix Z £ i^'^x" and t it holds that 

\\z\l<- 

II II* - 2 



if and only if 



3 symmetric matrices V G S™^™-,!^ G S"^" 
s.t. w)-^ "^'^ Tr(y) + TV(VF) = t 



Proof. =^ From Lemma 27 we obtain a matrix =:X^Oof trace say s < t. If s < t, 

we add (t — s) to the top-left entry of V , i.e. we add to X the PSD rank-1 matrix (t — s)eie^ 
(which again gives a PSD matrix). <^ follows directly from Lemma 27. □ 



11.2.1 Weighted Nuclear Norm 

A promising weighted nuclear norm regularization for matrix completion was recently proposed 
by [SSIO]. For fixed weight vectors p G M™,*/ G M", the weighted nuclear norm ||-Z^||n„c(p 9) 
^ e Rrnxn -g defined as 

where P = (i\&g{yfp) G M"^x™- denotes the diagonal matrix whose i-th diagonal entry is y/pl, 
and analogously for Q = diag(y^) G M"^". Here p G M"* is the vector whose entries are the 
probabilities p{i) > that the i-th row is observed in the sampling fi. Analogously, q G M" 
contains the probability q{j) > for each column j. The opposite weighting (using ^ and 
instead of p{i),q{j)) has also been suggested by [WKS08]. 

Any optimization problem with a weighted nuclear norm regularization 

min f(Z) (37) 

and arbitrary loss function / can therefore be formulated equivalently over the domain ||PZQ||^ < 
t/2, such that it reads as (if we substitute Z := PZQ), 

min f{p-^ZQ-^). 

ZeM™X",||z|| <t/2 

Hence, we have reduced the task to our standard convex problem (32) for / that here is defined 
as 

f{X) :=f{p-'ZQ-'), 

where X =: (^^). This equivalence implies that any algorithm solving (32) also serves as an 
algorithm for weighted nuclear norm regularization. In particular, Hazan's Algorithm 6 does 
imply a guaranteed approximation quality of e for problem (37) after O(^) many rank-1 updates, 
as we discussed in Section 7. So far, to the best of our knowledge, no approximation guarantees 
were known for the weighted nuclear norm. 

Solution path algorithms (maintaining approximation guarantees when the regularization 
parameter t changes) as proposed by [GJLIO, GJL12a, GJL12b], and the author's PhD the- 
sis [Jagll], can also be extended to the case of the weighted nuclear norm. 
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11.3 The Max-Norm for Matrices 



We think of the matrix max-norm as a generahzation of the vector ^oo-norm to the case of 
positive semidefinite matrices, which we have studied before. 

In some matrix completion apphcations, the max-norm has been observed to provide solutions 
of better generalization performance than the nuclear norm [SRJ04] . Both matrix norms can be 
seen as a convex surrogate of the rank [SS05]. 

The max-norm \\Z\\^^^ of a rectangular matrix Z G l^^x" has a nice characterization in terms 
of matrix factorizations of Z, i.e. 

Il^llmax := ^min^max{||L|||^ , \\R\\l^} , (38) 

where the number of columns of the factors L S MJ^^^ and R^^^ is not constrained [LRS^IO]. 
Here ||-^^||2oo maximum ^2-norm of any row Li- of L, that is ||-^^||2oo maxj ||Lj:||2 = 

\l^k ^"ik- Compared to the nuclear norm, we therefore observe that the max-norm con- 
strains the maximal Euclidean row-norms of any factorization of the original matrix Z, see 
also [SS05]. ^ 

An alternative formulation of the max-norm was given by [LMSS07] and [SS05], stating that 

||2'|| = min (max ||Lj:||2)(max | I2) . 

The dual norm to the max-norm, as given in [SS05], is 

ll-^llmax = Z»Y 

II Umax ii^ii 

I Umax — 

max ^ Zij ifrj , 
;,eM*,||ii||2<i 

r,eMM|r,||2<l 

where the last equality follows from the characterization (38). 

The relation of the max-norm to semidefinite optimization — which also explains the naming 
of the max-norm — is that 

ll-2^IL=^ = minimize t 
v,w 

. ^^ A Vie [m], (39) 

[Z^ w)"^^ Wu<t V^GN 

Here the two optimization variables range over the symmetric matrices V G g^xm 
W G S"^", see for example [LRS^IO]. As already in the nuclear norm case, this semidefi- 
nite characterization will again be the central tool for our algorithmic approach for max-norm 
regularized problems in the following. The equivalence of the above characterization to the 
earlier "factorization" formulation (38) is a consequence of the following simple Lemma 29. The 
Lemma gives a correspondence between the (rectangular) matrices Z £ ]R»^x" of bounded max- 
norm on one hand, and the (symmetric) PSD matrices X G g(m+n)x(m+n) q£ uniformly bounded 
diagonal on the other hand. 

Lemma 29. For any non-zero matrix Z € R^^m ^ ^ ]^_. 

ll^llmax<i 

^Note that the max-norm does not coincide with the matrix norm induced by the vector ||.||^-norm, that is 

II Zx I 

1 1 -^11 00 suPi/o ii^ii ■ The latter matrix norm by contrast is known to be the maximum of the row sums 
of Z (i.e. the £i-norms of the rows). 
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if and only if 



"nxn 



3 symmetric matrices V G g"ix»«^ G S 

V z\ , Vii<t Vie H, 



W J Wii<t [n] 



max 



Proof. =^ Using the above characterization (38) of the max-norm, or namely that \\Z\ 
minj^^T=2 max{||L||2 ^ , \\R\\2 ooi' '^^ S®* ^^^^ ^ with LR^ = Z, s.t. max{||L||2 ^ , II-RII2 ooi 
max{maxj ||-Z^i:||2 ,niaxj ||i?i;||2} ^ t, or in other words we have found a matrix {^^t /j^t ) = 
{L-R){L-RY h where every diagonal element is at most t, that is ||i^j:||2 = {LL^)ii <t\li£ 
[m], and \\Ri:\\l = {RR^h < i Vi G [n]. 

I I As the matrix ( ) is symmetric and PSD, it can be (Cholesky) factorized to (L; R){L; i?)^ 
s.t. Li?'^ = Z and ||i^i:||2 = {LL^)ii < t \/i £ [m] and ||-Ri:||2 = {RR^)ii <t\li£ [n], which 
implies ll^lUax < n 

11.4 Optimizing with Bounded Nuclear Norm and Max-Norm 

Most of the currently known algorithms for matrix factorizations as well as nuclear norm or 
max-norm regularized optimization problems, such as (31), (32), (33) or (34), do suffer from the 
following problem: 

In order to optimize the convex objective function f{Z) while controlling the norm or 
ll-^llmax' methods instead try to optimize f{LR^), with respect to both factors L G ig'^x's 
and R G M"'^'^, with the corresponding regularization constraint imposed on L and R. This 
approach is of course very tempting, as the constraints on the factors — which originate from 
the matrix factorization characterizations (35) and (38) — are simple and in some sense easier 
to enforce. 



Unhealthy Side-Effects of Factorizing. However, there is a significant price to pay: Even if the 
objective function f{Z) is convex in Z, the very same function expressed as a function f{LR?") of 
both the factor variables L and R becomes a severely non-convex problem, naturally consisting 
of a large number of saddle-points (consider for example just the smallest case L,R £ M}^^ 
together with the identity function f{Z) = Z G M). 

The majority of the currently existing methods such as for example [RS05, Lin07] and later also 
[Pat07, ZWSP08, KBV09, TPNT09, IRIO, GNHSll] is of this "alternating" gradient descent 
type, where at each step one of the factors L and R is kept fixed, and the other factor is 
updated by e.g. a gradient or stochastic gradient step. Therefore, despite working well in many 
practical applications, all these mentioned methods can get stuck in local minima — and so are 
theoretically not well justified, see also the discussion in [DeC06]. 

The same issue also comes up for max-norm optimization, where for example [LRS^IO] opti- 
mize over the non-convex factorization (38) for bounded max-norm. 

Concerning the fixed rank of the factorization, [GGIO] have shown that finding the optimum 
under a rank constraint (even if the rank is one) is NP-hard (here the used function / was the 
standard squared error on an incomplete matrix). On the positive side, [BM03] have shown 
that if the rank k of the factors L and R exceeds the rank of the optimum solution X* , then — 
in some cases — it can be guaranteed that the local minima (or saddle points) are also global 
minima. However, in nearly all practical applications it is computationally infeasible for the 
above mentioned methods to optimize with the rank k being in the same order of magnitude as 
the original matrix size m and n (as e.g. in the Netflix problem, such factors L, R could possibly 
not even be stored on a single machine*^). 

* Algorithm 6 in contrast does never need to store a full estimate matrix X, but instead just keeps the rank-1 
factors V obtained in each step, maintaining a factorized representation of X. 
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Relief: Optimizing Over an Equivalent Convex Problem. Here we simply overcome this 
problem by using the transformation to semidefinite matrices, which we have outlined in the 
above Corollary 28 and Lemma 29. These bijections of bounded nuclear and max-norm matrices 
to the PSD matrices over the corresponding natural convex domains do allow us to directly 
optimize a convex problem, avoiding the factorization problems explained above. We describe 
this simple trick formally in the next two Subsections 11.4.1 and 11.4.2. 



But what if you really need a Matrix Factorization? In some applications (such as for exam- 
ple embeddings or certain collaborative filtering problems) of the above mentioned regularized 
optimization problems over f{Z), one would still want to obtain the solution (or approximation) 
Z in a factorized representation, that is Z = LR^ . We note that this is also straight-forward to 
achieve when using our transformation: An explicit factorization of any feasible solution to the 
transformed problem (20) or (27) — if needed — can always be directly obtained since X ^ 0. 

Alternatively, algorithms for solving the transformed problem (20) can directly maintain the 
approximate solution X in a factorized representation (as a sum of rank-1 matrices), as achieved 
for example by Algorithms 6 and 7. 



11.4.1 Optimization with a Nuclear Norm Regularization 

Having Lemma 27 at hand, we immediately get to the crucial observation of this section, allowing 
us to apply Algorithm 6: 

Any optimization problem over bounded nuclear norm matrices (32) is in fact equivalent to 
a standard bounded trace semidefinite problem (20). The same transformation also holds for 
problems with a bound on the weighted nuclear norm, as given in (37). 

Corollary 30. Any nuclear norm regularized problem of the form (32) is equivalent to a bounded 
trace convex problem of the form (20), namely 

minimize f{X) 

s.t. Tr{X) = t , (40) 
X hO 

where f is defined by f{X) := f{Z) for Z G K™x" being the upper right part of the symmetric 
matrix X. Formally we again think of X G g(n+m)x(n+m) consisting of the four parts X =: 

xf. ) with V G S'^xm^ ^ g gnxn ^ ^ M'"Xn_ 

Here "equivalent" means that for any feasible point of one problem, we have a feasible point 
of the other problem, attaining the same objective value. The only difference to the original 
formulation (20) is that the function argument X needs to be rescaled by | in order to have 
unit trace, which however is a very simple operation in practical applications. Therefore, we can 
directly apply Hazan's Algorithm 6 for any max-norm regularized problem as follows: 

Using our analysis of Algorithm 6 from Section 7.1, we see that Algorithm 8 runs in time near 
linear in the number iVj of non-zero entries of the gradient V/. This makes it very attractive in 
particular for recommender systems applications and matrix completion, where V/ is a sparse 
matrix (same sparsity pattern as the observed entries), which we will discuss in more detail in 
Section 11.5. 

Corollary 31. After at most O (^) many iterations (i.e. approximate eigenvalue computations), 
Algorithm 8 obtains a solution that is e close to the optimum of (32). The algorithm requires a 
total of O arithmetic operations (with high probability). 
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Algorithm 8 Nuclear Norm Regularized Solver 

Input: A convex nuclear norm regularized problem (32), 

target accuracy e 
Output: e-approximate solution for problem (32) 

1. Consider the transformed symmetric problem for /, 

as given by Corollary 30 

2. Adjust the function / so that it first rescales its argument by t 

3. Run Kazan's Algorithm 6 for f{X) over the domain X £ S. 



Proof. We use the transformation from Corollary 30 and then rescale all matrix entries by 
J. Then result then follows from Corollary 20 on page 31 on the running time of Kazan's 
algorithm. □ 

The fact that each iteration of our algorithm is computationally very cheap — consisting only 
of the computation of an approximate eigenvector — strongly contrasts the existing "proximal 
gradient" and "singular value thresholding" methods [GLW+09, JY09, MGC09, LST09, CCSIO, 
TYIO], which in each step need to compute an entire SVD. Such a single incomplete SVD 
computation (first k singular vectors) amounts to the same computational cost as an entire 
run of our algorithm (for k steps). Furthermore, those existing methods which come with a 
theoretical guarantee, in their analysis assume that all SVDs used during the algorithm are 
exact, which is not feasible in practice. By contrast, our analysis is rigorous even if the used 
eigenvectors are only e'- approximate. 

Another nice property of Kazan's method is that the returned solution is guaranteed to be 
simultaneously of low rank (k after k steps), and that by incrementally adding the rank-1 
matrices Vkv'^ , the algorithm automatically maintains a matrix factorization of the approximate 
solution. 

Also, Kazan's algorithm, as being an instance of our presented general framework, is designed 
to automatically stay within the feasible region S, where most of the existing methods do need 
a projection step to get back to the feasible region (as e.g. [Lin07, LST09]), making both the 
theoretical analysis and implementation more complicated. 

11.4.2 Optimization with a Max-Norm Regularization 

The same approach works analogously for the max-norm, by using Lemma 29 in order to apply 
Algorithm 7: 

Any optimization problem over bounded max-norm matrices (34) is in fact equivalent to a 
semidefinite problem (27) over the "box" of matrices where each element on the diagonal is 
bounded above by t. We think of this domain as generalizing the positive cube of vectors, to 
the PSD matrices. 

Corollary 32. Any max-norm regularized problem of the form (34) is equivalent to a bounded 
diagonal convex problem of the form (27), i.e., 

minimize fiX) 

XG§(™+") X (.rn+n) 

s.t. Xii < 1 yi, (41) 
XhO 

where f is defined by f{X) := f{Z) for Z G i^™x" heing the upper right part of the symmetric 
matrix X. Formally we again think of any X S g(n+m)x(n+m) consisting of the four parts 

X =: (^j, ^ j with V £ S™^"^, W £ S"^" and Z £ M'"^". 
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Again the only difference to the original formulation (27) is that the function argument X 
needs to be rescaled by j in order to have the diagonal bounded by one, which however is a 
very simple operation in practical applications. This means we can directly apply Algorithm 7 
for any max-norm regularized problem as follows: 

Algorithm 9 Max-Norm Regularized Solver 

Input: A convex max-norm regularized problem (34), 

target accuracy e 
Output: e-approximate solution for problem (34) 

1. Consider the transformed symmetric problem for /, 

as given by Corollary 32 

2. Adjust the function / so that it first rescales its argument by t 

3. Run Algorithm 7 for f{X) over the domain X G ffl. 



Using the analysis of our new Algorithm 7 from Section 7.1, we obtain the following guarantee: 



8Cf 



many iterations, Algorithm 9 obtains a solution that is £ close to the 



Corollary 33. After 
optimum of (34)- 

Proof. We use the transformation from Corollary 32 and then rescale all matrix entries by j. 
Then the running time of the algorithm follows from Theorem 25. □ 



Maximum Margin Matrix Factorizations. In the case of matrix completion, the "loss" func- 
tion / is defined as measuring the error from X to some fixed observed matrix, but just at a 
small fixed set of "observed" positions of the matrices. As we already mentioned, semidefinite 
optimization over X as above can always be interpreted as finding a matrix factorization, as a 
symmetric PSD matrix X always has a (unique) Cholesky factorization. 

Now for the setting of matrix completion, it is known that the above described optimization 
task under bounded max-norm, can be geometrically interpreted as learning a maximum margin 
separating hyperplane for each user/movie. In other words the factorization problem decomposes 
into a collection of SVMs, one for each user or movie, if we think of the corresponding other 
factor to be fixed for a moment [SRJ04]. We will discuss matrix completion in more detail in 
Section 11.5. 



Other Applications of Max-Norm Optimization. Apart from matrix completion, optimiza- 
tion problems employing the max-norm have other prominent applications in spectral methods, 
spectral graph properties, low-rank recovery, and combinatorial problems such as Max-Cut. 



11.5 Applications 

Our Algorithm 8 directly applies to arbitrary nuclear norm regularized problems of the form (32). 
Since the nuclear norm is in a sense the most natural generalization of the sparsity-inducing ii- 
norm to the case of low rank matrices (see also the discussion in the previous chapters) there 
are many applications of this class of optimization problems. 



11.5.1 Robust Principal Component Analysis 

One prominent example of a nuclear norm regularized problem in the area of dimensionality 
reduction is given by the technique of robust PC A as introduced by [CLMWll], also called 
principal component pursuit, which is the optimization task 

min ||Z||^ +/x||M-Z||^ . (42) 



m X n 
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Here M £ j^'^x'^ the given data matrix, and ||.||-|^ denotes the entry- wise £i-norm. By consid- 
ering the equivalent constrained variant \\Z\\^ < ^ instead, we obtain a problem the form (32), 
suitable for our Algorithm 8. However, since the original objective function f{Z) = \\M — Z\\-^ 
is not differentiable, a smoothed version of the £i-norm has to be used instead. This situation 
is analogous to the hinge-loss objective in maximum margin matrix factorization [SRJ04]. 

Existing algorithms for robust PCA do usually require a complete (and exact) SVD in each 
iteration, as e.g. [TYIO, AGIll], and are often harder to analyze compared to our approach. The 
first algorithm with a convergence guarantee of O(^) was given by [AGIll], requiring a SVD 
computation per step. Our Algorithm 8 obtains the same guarantee in the same order of steps, 
but only requires a single approximate eigenvector computation per step, which is significantly 
cheaper. 

Last but not least, the fact that our algorithm delivers approximate solutions to (42) of rank 
O (^) will be interesting for practical dimensionality reduction applications, as it re-introduces 
the important concept of low-rank factorizations as in classical PCA. In other words our algo- 
rithm produces an embedding into at most O (^) many new dimensions, which is much easier 
to deal with in practice compared to the full rank n solutions resulting from the existing solvers 
for robust PCA, see e.g. [CLMWll] and the references therein. 

We did not yet perform practical experiments for robust PCA, but chose to demonstrate the 
practical performance of Algorithm 6 for matrix completion problems first. 

11.5.2 Matrix Completion and Low Norm Matrix Factorizations 

For matrix completion problems as for example in collaborative filtering and recommender sys- 
tems [KBV09], our algorithm is particularly suitable as it retains the sparsity of the observations, 
and constructs the solution in a factorized way. In the setting of a partially observed matrix 
such as in the Netflix case, the loss function f{X) only depends on the observed positions, which 
are very sparse, so Vf{X) — which is all we need for our algorithm — is also sparse. 

We want to approximate a partially given matrix Y (let be the set of known training entries 
of the matrix) by a product Z = LR^ such that some convex loss function f{Z) is minimized. 
By iltest we denote the unknown test entries of the matrix we want to predict. 

Complexity. Just recently it has been shown that the standard low-rank matrix completion 
problem — that is finding the best approximation to an incomplete matrix by the standard £2- 
norm — is an NP-hard problem, if the rank of the approximation is constrained. The hardness 
is claimed to hold even for the rank 1 case [GGIO]. In the light of this hardness result, the 
advantage of relaxing the rank by replacing it by the nuclear norm (or max-norm) is even 
more evident. Our near linear time Algorithm 8 relies on a convex optimization formulation 
and does indeed deliver an guaranteed e-accurate solution for the nuclear norm regularization, 
for arbitrary e > 0. Such a guarantee is lacking for the "alternating" descent heuristics such 
as [RS05, Lin07, Pat07, ZWSP08, KBV09, TPNT09, IRIO, GNHSll, SSIO, LRS+10, RRll] 
(which build upon the non-convex factorized versions (35) and (38) while constraining the rank 
of the used factors L and R) . 

Different Regularizations. Regularization by the weighted nuclear norm is observed by [SSIO] 
to provide better generalization performance than the classical nuclear norm. As it can be simply 
reduced to the nuclear norm, see Section 11.2.1, our Algorithm 8 can directly be applied in the 
weighted case as well. On the other hand, experimental evidence also shows that the max-norm 
sometimes provides better generalization performance than the nuclear norm [SRJ04, LRS+10]. 
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For any convex loss function, our Algorithm 9 solves the corresponding max-norm regularized 
matrix completion task. 

Different Loss Functions. Our method applies to any convex loss function on a low norm 
matrix factorization problem, and we will only mention two loss functions in particular: 

Maximum Margin Matrix Factorization (MMMF) [SRJ04] can directly be solved by our Al- 
gorithm 8. Here the original (soft margin) formulation is the trade-off formulation (31) with 
f{Z) := X^jjgQ \ Zij — yij\ being the hinge or £i-loss. Because this function is not differentiable, 
the authors recommend using the differentiable smoothed hinge loss instead. 

When using the standard squared loss function f{Z) := J2ijen(^v ~ Vv)"^ ^ problem is 
known as Regularized Matrix Factorization [Wu07], and both our algorithms directly apply. 
This loss function is widely used in practice, has a very simple gradient, and is the natural 
matrix generalization of the ^2-loss (notice the analogous Lasso and regularized least squares 
formulation). The same function is known as the rooted mean squared error, which was the 
quality measure used in the Netflix competition. We write RMSEtrain and RMSEfest for the 
rooted error on the training ratings and test ratings litest respectively. 

Running time and memory. From Corollary 31 we have that the running time of our nuclear 
norm optimization Algorithm 8 is linear in the size of the input: Each matrix- vector multipli- 
cation in Lanczos' or the power method exactly costs \^}\ (the number of observed positions 
of the matrix) operations, and we know that in total we need at most O (l/e^'^) many such 
matrix- vector multiplications. 

Also the memory requirements are very small: Either we store the entire factorization of X^'^'^ 
(meaning the O (-) many vectors v^''^) — which is still much smaller than the full matrix X 

(k) 

— or then instead we can only update and store the prediction values X^ for ij G 17 U fltest 
in each step. This, together with the known ratings yij determines the sparse gradient matrix 
Wf{X^^^) during the algorithm. Therefore, the total memory requirement is only |r2ur2iest| 
(the size of the output) plus the size (n -|- m) of a single feature vector v. 

The constant in the running time of Algorithm 6. One might ask if the constant hidden 
in the O(^) number of iterations is indeed controllable. Here we show that for the standard 
squared error on any fixed set of observed entries this is indeed the case. For more details on 
the constant C/, we refer the reader to Sections 3.4 and 7.1. 

Lem ma 34. For the squared error f{Z) — ^ ^^/ij^fii^ij — over the spectahedron S, it holds 
that Cj <1. 

Proof. In Lemma 6, we have seen that the constant Cj is upper bounded by half the diameter 
of the domain, times the largest eigenvalue of the Hessian V"^ f{X). Here we consider / as a 

— * 2 

function on vectors X S corresponding to the matrices X G S"^". However for the squared 
error as in our case here, the Hessian will be a diagonal matrix. One can directly compute that 
the diagonal entries of V^/(X) are 1 at the entries corresponding to 12, and zero everywhere 
else. Furthermore, the squared diameter of the spectahedron is upper bounded by 2, as we have 
shown in Lemma 23. Therefore < 1 for the domain S. □ 

If the domain is the scaled spectahedron t • 5 as used in our Algorithm 8, then the squared 
diameter of the domain is 2t^, compare to Lemma 23. This means that the curvature is upper 
bounded by < in this case. Alternatively, the same bound for the curvature of f{X) := 

f{tX) can be obtained along the same lines as for the spectahedron domain in the previous 
lemma, and the same factor of t"^ will be the scaling factor of the Hessian, resulting from the 
chain-rule for taking derivatives. 
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11.5.3 The Structure of the Resulting Eigenvalue Problems 

For the actual computation of the approximate largest eigenvector in Algorithm 6, i.e. the 

— V/(X('^')), jrp^ j , either Lanczos' method or the power method 
(as in PageRank, see e.g. [Ber05]) can be used. In our Theorem 18 of Section 7.1, we stated 
that both the power method as well as Lanczos' algorithm do provably obtain the required 
approximation quality in a bounded number of steps if the matrix is PSD, with high probability, 
see also [KW92, AHK05]. 

Both methods are known to scale well to very large problems and can be parallelized easily, 
as each iteration consists of just one matrix-vector multiplication. However, we have to be 
careful that we obtain the eigenvector for the largest eigenvalue which is not necessarily the 
principal one (largest in absolute value). In that case the spectrum can be shifted by adding an 
appropriate constant to the diagonal of the matrix. 

For arbitrary loss function /, the gradient — V/(X), which is the matrix whose largest eigen- 
vector we have to compute in the algorithm, is always a symmetric matrix of the block form 

Vf{X) = for G = Vf{Z), when X =: In other words V/(X) is the 

adjacency matrix of a weighted bipartite graph. One vertex class corresponds to the n rows of 
the original matrix X2 {users in recommender systems), the other class corresponds to the m 
columns {items or movies). It is easy to see that the spectrum of V/ is always symmetric: 
Whenever (^) is an eigenvector for some eigenvalue A, then {-w) is an eigenvector for —A. 

Hence, we have exactly the same setting as in the established Hubs and Authorities (HITS) 
model [Kle99]. The first part of any eigenvector is always an eigenvector of the hub matrix G^G, 
and the second part is an eigenvector of the authority matrix GG^ . 

Repeated squaring. In the special case that the matrix G is very rectangular (n ^ m or 
n ^ m), one of the two square matrices G^G or GG^ is very small. Then it is known that one 
can obtain an exponential speed-up in the power method by repeatedly squaring the smaller 
one of the matrices, analogously to the "square and multiply" -approach for computing large 
integer powers of real numbers. In other words we can perform 0(log ^) many matrix-matrix 
multiplications instead of O(^) matrix-vector multiplications. 

11.5.4 Relation to Simon Funk's SVD Method 

Interestingly, our proposed framework can also be seen as a theoretically justified variant of 
Simon Funk's [Wcb06] and related approximate SVD methods, which were used as a building 
block by most of the teams participating in the Netflix competition (including the winner team). 
Those methods have been further investigated by [Pat07, TPNT09] and also [KBC07], which 
already proposed a heuristic using the HITS formulation. These approaches are algorithmically 
extremely similar to our method, although they are aimed at a slightly different optimization 
problem, and do not directly guarantee bounded nuclear norm. Recently, [SSIO] observed that 
Funk's algorithm can be seen as stochastic gradient descent to optimize (31) when the regular- 
ization term is replaced by a weighted variant of the nuclear norm. 

Simon Funk's method considers the standard squared loss function f{X) = | '^ij^zsi-^ii ~ 
Uij)'^, and finds the new rank-1 estimate (or feature) v by iterating v := v + X{—Vf{X)v — Kv), 
or equivalently 

v:=X (-V/(X) +(^^-K^I^v, (43) 

a fixed number of times. Here A is a small fixed constant called the learning rate. Additionally 
a decay rate K > is used for regularization, i.e. to penalize the magnitude of the resulting 
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feature v. This matrix- vector multiplication formulation (43) is equivalent to a step of the power 
method applied within our framework^, and for small enough learning rates A the resulting 
feature vector will converge to the largest eigenvector of — V/(Z). 

However in Funk's method, the magnitude of each new feature strongly depends on the starting 
vector vq, the number of iterations, the learning rate A as well as the decay K, making the 
convergence very sensitive to these parameters. This might be one of the reasons that so far no 
results on the convergence speed could be obtained. Our method is free of these parameters, the 
k-th new feature vector is always a unit vector scaled by Also, we keep the Frobenius norm 

1 1 2 1 1 2 

ll^llFro + ll^llFro obtained factorization exactly fixed during the algorithm, whereas in 

Funk's method — which has a different optimization objective — this norm strictly increases 
with every newly added feature. 

Our described framework therefore gives theoretically justified variant of the experimentally 
successful method [Wcb06] and its related variants such as [KBC07, Pat07, TPNT09]. 

11.6 Experimental Results 

We run our algorithm for the following standard datasets^*^ for matrix completion problems, 
using the squared error function. 



dataset 


#ratings 


n 


m 


MovieLens 100k 


10^ 


943 


1682 


MovieLens IM 


10*^ 


6040 


3706 


MovieLens lOM 


10^ 


69878 


10677 


Netflix 




480189 


17770 



Any eigenvector method can be used as a black-box in our algorithm. To keep the experiments 
simple, we used the power method^^, and performed 0.2 • k power iterations in step k. If not 
stated otherwise, the only optimization we used is the improvement by averaging the old and 
new gradient as explained in Section 7.3. All results were obtained by our (single-thread) 
implementation in Java 6 on a 2.4 GHz Intel C2D laptop. 

Sensitivity. The generalization performance of our method is relatively stable under changes 
of the regularization parameter, see Figure 2: 

Movielens. Table 1 reports the running times of our algorithm on the three MovieLens datasets. 
Our algorithm gives an about 5.6 fold speed increase over the reported timings by [TYIO], which 
is a very similar method to [JY09]. [TYIO] already improves the "singular value thresholding" 
methods [CCSIO] and [MGC09]. For MMMF, [RS05] report an optimization time of about 5 
hours on the IM dataset, but use the different smoothed hinge loss function so that the results 
cannot be directly compared. [MGC09], [SJ03] and [JY09] only obtained results on much smaller 
datasets. 

In the following experiments on the MovieLens and Netflix datasets we have pre- normalized 
all training ratings to the simple average ^' of the user and movie mean values, for the sake 
of being consistent with comparable literature. 

^Another difference of our method to Simon Funk's lies in the stochastic gradient descent type of the latter, 
i.e. "immediate feedback" : During each matrix multiplication, it already takes the modified current feature 
V into account when calculating the loss f{Z), whereas our Algorithm 6 alters Z only after the eigenvector 
computation is finished. 
^''See www.grouplcns.org and archive. ics.uci.cdu/ml. 

""^^ We used the power method starting with the uniform unit vector. | of the approximate eigenvalue corresponding 
to the previously obtained feature Vk-i was added to the matrix diagonal to ensure good convergence. 
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0.95 



0.93 



0.91 




0.89 

15000 30000 45000 60000 

Trace regularization t 

Figure 2: Sensitivity of the method on the choice of the regularization parameter t in (32), on MovieLens 
IM. 

Table 1: Running times iour (in seconds) of our algorithm on the three MovieLens datasets compared to 
the reported timings ^ty of [TYIO]. The ratings {1, . . . , 5} were used as-is and not normalized 
to any user and/or movie means. In accordance with [TYIO], 50% of the ratings were used 
for training, the others were used as the test set. Here NMAE is the mean absolute error, 
times 5ir[, over the total set of ratings, k is the number of iterations of our algorithm, ^mm 
is the total number of sparse matrix-vector multiplications performed, and tr is the used trace 
parameter t in (32). They used Matlab/PROPACK on an Intel Xeon 3.20 GHz processor. 





NMAE 


ixY 


^our 


k 




tr 


100k 


0.205 


7.39 


0.156 


15 


33 


9975 


IM 


0.176 


24.5 


1.376 


35 


147 


36060 


lOM 


0.164 


202 


36.10 


65 


468 


281942 



For MovieLens lOM, we used partition r^, provided with the dataset (10 test ratings per user). 
The regularization parameter t was set to 48333. We obtained a RMSEtest of 0.8617 after k = 400 
steps, in a total running time of 52 minutes (16291 matrix multiplications). Our best RMSEtest 
value was 0.8573, compared to 0.8543 obtained by [LU09] using their non-linear improvement 
of MMMF. 

Algorithm Variants. Comparing the proposed algorithm variants from Section 7.3, Figure 3 
demonstrates moderate improvements compared to our original Algorithm 8. 

Netflix. Table 2 compares our method to the two "hard impute" and "soft impute" sin- 
gular value thresholding methods of [MHTIO] on the Netflix dataset, where they used Mat- 
lab/PROPACK on an Intel Xeon 3 GHz processor. The "soft impute" variant uses a constrained 
rank heuristic in each update step, and an "un-shrinking" or fitting heuristic as post-processing. 
Both are advantages for their method, and were not used for our implementation. Nevertheless, 
our algorithm seems to perform competitive compared to the reported timings of [MHTIO]. 

Note that the primary goal of this experimental section is not to compete with the prediction 
quality of the best engineered recommender systems (which are usually ensemble methods, i.e. 
combinations of many different individual methods). We just demonstrate that our method 
solves nuclear norm regularized problems of the form (32) on large sample datasets, obtaining 
strong performance improvements. 
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0.94 



0.863 



0.785 



RMSE 



0.708 




0.63 



100 



200 
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Figure 3: Improvements for the two algorithm variants described in Section 7.3, when running on Movie- 
Lens lOM. The thick hnes above indicate the error on the test set, while the thinner lines 
indicate the training error. 



Table 2: Running times tour (in hours) of our algorithm on the Netflix dataset compared to the reported 
timings iivi.hard for "hard impute" by [MHT09] and tM,soft for "soft impute" by [MHTIO]. 



RMSEtest 


tM,hard 


tlVljSoft 


tour 


k 


#mm 


tr 


0.986 


3.3 


n.a. 


0.144 


20 


50 


99592 


0.977 


5.8 


n.a. 


0.306 


30 


109 


99592 


0.965 


6.6 


n.a. 


0.504 


40 


185 


99592 


0.962 


n.a. 


1.36 


1.08 


45 


243 


174285 


0.957 


n.a. 


2.21 


1.69 


60 


416 


174285 


0.954 


n.a. 


2.83 


2.68 


80 


715 


174285 


0.9497 


n.a. 


3.27 


6.73 


135 


1942 


174285 


0.9478 


n.a. 


n.a. 


13.6 


200 


4165 


174285 



11.7 Conclusion 



We have introduced a new method to solve arbitrary convex problems with a nuclear norm 
regularization, which is simple to implement and to parallelize. The method is parameter- 
free and comes with a convergence guarantee. This guarantee is, to our knowledge, the first 
guaranteed convergence result for the class of Simon- Funk-type algorithms, as well as the first 
algorithm with a guarantee for max-norm regularized problems. 

It remains to investigate if our algorithm can be applied to other matrix factorization problems 
such as (potentially only partially observed) low rank approximations to kernel matrices as used 
e.g. by the PSVM technique [CZW^OT], regularized versions of latent semantic analysis (LSA), 
or non-negative matrix factorization [Wu07]. 
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